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ABSTRACT 


In this thesis, a fast algorithm for the parameter 
identification of linear stochastic systems is derived and 


its unbiased convergence property is proved. 


To lay the ground work for the above development, a 
new faster algorithm for linear deterministic System is 
derived by improving the previously published Kudva-Narendra 
algorithm. Geometric interpretation of the Kudva-Narendra 
algorithm helps one to understand how this algorithm works 
and consequently to improve the convergence rate, thus 
resulting in a rapidly converging algorithm with reduced 


computational requirements. 


The algorithm for the stochastic system 1s obtained by 
modifying the Extended Kalman Filter algorithm using a 
canonical form of the innovations representation and using 
the concepts developed for the deterministic case. This 
development reveals that other existing algorithms which 
were developed using different approaches can be shown to be 


modifications of the EKF algorithm. 


The effectiveness of this algorithm is compared with 
two other existing algorithms by means of simulation 
Studies. These studies show the superiority of the new 
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algorithm with respect to convergence rate and computational 


requirements. 
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CHAPTER I 
GENERAL INTRODUCTION 


lm introduction 

Any successful application of modern control theory 
largely depends on the information available with regard to 
the physical system of interest. This information is 
usually incorporated into a mathematical model that behaves 
identical to the physical system in terms of its 
input-output relations. The mathematical model can be 
characterized mainly by two aspects: 

as istructure 

b. parameters 
Therefore, the collection of the information about the 
physical system is equivalent to the determination of the 
structure and parameters of the mathematical model on the 
basis of input and output information (Zadeh, 1962). Sucha 
determination is often referred to as the identification 


problem and is the main subject pursued in this thesis. 


It is extremely useful if it is possible to identify a 
black box only with input and output information. However, 
a baSic principle of identification is that a completely 
unknown system cannot be identified. The first step, then, 
to be taken in the formulation of an identification problem 
is to characterize the system. This characterization, in 


general, involves the selection of an appropriate structure 
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for the system.,” The structure chosen should reflect factors 
such as, any a priori knowledge of the system and the 


ultimate purpose of the identification. 


Due. to the advent of the digital computer, it is now 
possible to employ sophisticated methods to control systems. 
Such methods require, however, that on-line estimation of 
the parameters be made even as the implementation of the 
control scheme is in progress. Recursive estimation is 
often mandatory in many modern adaptive control systems. 
Thus, considerable attention has recently been paid to 


recurSive identification algorithms. 


Of the many identification techniques which have been 
developed in recent years, the model reference technique 
possesses a strong intuitive appeal in that it fully 
utilizes the advantage of the recursive algorithm. The use 
of the term model reference in the context of identification 
perhaps needs some justification. As first conceived by 
Whittaker (1958) for aircraft control, the model reference 
technique adjusts the controller parameters in such a manner 
that the parameters of the composite system namely, the 
controller and the system to be controlled, approach the 
parameters of the reference model in some optimal sense. 

The philosophy of this technique can be used for the 
identification problem by adjusting the model parameters to 
follow the parameters of the unknown system, which is looked 


upon as the reference system. This adjustment is best done 
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inva recursive fashion. A natural way to come up with a 
Suitable recursive algorithm would be to introduce a cost 
function in terms of the error between the mode? and the 
system, and to devise an adjustment mechanism for minimizing 
the-cost function. This approach not only converts the 
problem of the model reference identification into an 
optimization problem but also provides a quantitative 


measure of the performance of the adjustment mechanism. 


The recurSive least-squares (LS) algorithm is perhaps 
the best known and widely used for linear systems (Astrom, 
et al., 1971). It is well known that this algorithm leads 
to a biased estimate when noise, if present, is correlated 
as 1s the case in the identification problem of linear 
stochastic systems. In order to avoid this serious drawback 
a number of modifications has been proposed (Isermann, 


“et -al.Rv1974; tSaridisy 1974 )Soderstrom, (etral.., 31978). 


Among them the Extended Least Squares (ELS) algorithm 
(Panuska, 1968: 1969) is the easiest to implement and 
consequently has gained popularity in real-life applications 
(Ljung, et al., 1975). The basic idea of this approach is 
that noise is regarded as an additional input signal and the 
LS technique is applied to it. Since the noise is not 
available for measurement, it is replaced by previous 
residuals. Early simulation studies seemed to indicate that 
this method had good convergence property (Panuska, 1968; 
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(1975: 1977), however, hadS revealed that, contrary to prior 
Simulation results, the ELS method will not always converge 
to the true values of the parameters suing ebony eecence is 
dependent on a condition predetermined by the unknown 
process parameters. This nonconvergence, which apparently 
stems from the replacement of noise by previous residuals, 
Gan) De corrected by the introduction of ay prefilter 
involving the unknown parameters. Approximating the 
prefilter recursively with the updated estimates results in 
the recurSive approximate maximum likelihood (RML) method 
(Soderstrom, 1973). This algorithm was developed by an 
approximation of the off-line maximum likelihood (ML) method 
(Astrom, et al., 1966) and it possesses a nice local 
convergence property. No global convergence results are, 


however, available at this time (Panuska, 1980b). 


By the inclusion of the unknown parameters into its 
state vector a linear stochastic system becomes a nonlinear 
Stochastic system. Hence the identification problem of a 
linear stochastic system is equivalent to the state 
estimation problem of a nonlinear stochastic system. Among 
the numerous approaches available for the nonlinear 
estimation problem, the Extended Kalman Filter (EKF) 
approach turns out to be suitable in many cases (Jazwinski, 
1970). Since the EKF approach is based on linearization at 
each step of iteration around the previous estimates and the 


application of the Kalman filtering theory (Kalman, 1960), 
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the algorithm usually suffers from the lack of global 
Stability and often involves excessive computation. This 
provides the motivation for modifying the EKF algorithm into 
purely a parameter identifier in order to reduce the 


computational requirements. 


The main difficulty associated with the stochastic 
identification problem has peer largely the lack of methods 
for the analysis of convergence properties. This motivated 
the extension of the stability theories, which had been 
developed mainly for the deterministic problem, to the 
Stochastic problem (Kushner, 1967; Jumarie, 1979). Of a 
number of stability theories, the hyperstability theory 
(Popov, 1962; Zames, 1966) has been most extensively used in 
various applications to the deterministic problem in order 
to obtain adaptive stable systems (Landau, 1974). Thus, it 
is a rather natural choice to extend the applicable area of 
the hyperstability theorem to the stochastic identification 
problem. This was first successfully achieved by Landau 
(1976) with an introduction of a decreaSing gain matrix to 
estimate the parameters of a linear system contaminated by a 
colored noise. This work was again extended by Ljung 
(1977a) to estimate the parameters of the noise as well as 
the system. The real importance of the Ljung's work, 
however, rests on finding the role of a positive real 
transfer function in the ELS method. This finding not only 


enables one to define the region of convergence but also 
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Suggests a way leading to a convergent algorithm. 


1.2 Scope of the Thesis 
The main objective of this thesis is the discussion 
and the development of several model reference adaptive 


identification algorithms for linear stochastic systems. 


In Chapter II, a representation of the linear 
Stochastic System 1s fully discussed in order to facilitate 
the formulation of the identification problem treated in 
Chapter IV. Possible approaches to be taken, in general, to 


this problem are also discussed. 


Chapter III presents a fast adaptive identification 
method for a linear deterministic system. The main reason 
for doing this is to develop certain ideas which will be 
used to solve the identification problem for the linear 
Stochastic system. The identification problem for the 
stochastic case is, in general, a nonlinear problem and not 
amenable to direct solution. Therefore, a solution will be 
attempted by linearization. In this task, the concepts used 
in the development of the algorithm for a deterministic 


System will be useful. 


In Section 3.2, a geometric interpretation of the 
Kudva-Narendra scheme (Kudva, et al., 1974) is given. This 
approach, which is very crucial to the development of the 
fast algorithm, provides greater insight into and 


consequently a clearer understanding of the K-N scheme. In 
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Secriones, 37 therconceptrotsorthogonaleqradienterspusedcin 
order to develop the algorithm for rapid convergence. In 
Section 3.4, the effectiveness of the algorithm proposed is 


demonstrated by the simulation studies 


Chapter IV deals with the linear stochastic 
identification problem. Since this identification problem 
as formulated is inherently a nonlinear estimation problem, 
approximation techniques are used to eorain a filter that is 
physically realizable. Among the many techniques that can 
be used, the EKF technique is of particular interest because 
Cleres Simplicity end thus populanmity in ‘real vappiications:. 
Section 4.2 discusses the application of the EKF to the 
identification problem. In an attempt to reduce the 
computational requirements which arise in this approach, 
modifications of the EKF algorithm are made resulting in two 
Simplified algorithms. These are given in Section 4.3 and 
4.4, For obvious reasons these are labelled as EKF-M1 and 
EKF-M2, respectively. Though these two algorithms have 
Similarities to some existing ones to the author's best 
knowledge, the exact algorithms reported here have not 


appeared in the literature. 


Convergence analysis is probably one of the most 
important means to assess the usefulness of identification 
algorithms. However, since EKF-M1 turns out to be similar 
to the RML algorithm, no detailed discussion of its 


convergence property is included here. In Section 4.5, 
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attention is therefore focused on the convergence property 
of EKF-M2. In Section 4.6, a generalization of EKF-M2 
leading to a third algorithm EKF-M3 is discussed. 

Comparison of the features of EKF-M3 with those of two other 
existing algorithms, namely, the RML method (Soderstrom, 


1973) and the ELS method (Panuska, 1968; 1969). 


The usefulness of EKF-M3 is evaluated by means of 
Simulation studies. For purpose of comparison, simulation 
Studies are also carried out with the RML and ELS 


algorithms. These resuls are reported in Chapter V. 


Chapter VI provides a summary of the thesis, 


conclusions and suggestions for further research. 
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CHAPTERS I 
FORMULATION OF IDENTIFICATION PROBLEM 


Zoleluncroduction 

An identification problem is characterized by three 
quantities (Zadeh, 1962): a class of model, a class of input 
Signals and a criterion. The selection of the class of 
models generally leads to the determination of the structure 
that mathematically represents the system. Since input 
Signals have a profound effect on the performance of the 
identification scheme some limitations are usually imposed 
on input signals. It has been shown (Aoki, et al., 1970; 
Yuan, et al., 1977) that the persistent excitation condition 
of input signals is sufficient in many cases to obtain 
consistent estimates or improve the performance of the 
identification scheme. However, the selection of input 
Signals itself is the subject of an interesting but separate 
area of research, hence we will not pursue this matter any 


further: 


The structure of models describing the system under 
consideration wields considerable influence on the 
complexity of identification schemes in various ways such as 
computational requirements, consistency and uniqueness of 
estimates, etc. Since different identification methods are 
necessary for different representations there are very few 


general rules available by which the optimality of the model 
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structure chosen can be determined. A corvenient 
characterization of the model for the recursive model 
reference approach to the identification problem can be made 
through the parametric representation since this is 


relatively easy to implement on computers. 


In the parametric representation of linear systems, 
the choice of model structure by which the system is 
represented is influenced by the nOmDer of unknown 
parameters to be identified and the applicability of the 
real-time optimization techniques which are readily 
available now. The number of parameters to be identified is 
directly related to the efficiency of the identification 
algorithm as well as the uniqueness of the solution because 
there exists a limit to number of the paramters which can be 
uniquely determined by the input-output relations. Since 
the identification problem is often formulated as an 
optimization problem by introducing a cost function, the 
applicability of the available optimization techniques 
should be taken into account in the choice of a model 


Structure. 


In this chapter we restrict our discussion to linear 
multiinput singleoutput (MISO) stochastic systems for the 
sake of simplicity. Since a linear multiinput multioutput 
(MIMO) system can be regarded as a group of MISO systems 
each of which has an input-output pair associated with it, 


the extension to MIMO systems can be readily carried out. 
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2.2 Representation of MISO Linear Stochastic Systems 

Kalman (1963) has shown that only the completely 
controllable and completely observable part of a system can 
be characterized from input-output relations. Thus, the 
system that 1s of interest in this thesis is one which is 
completely controllable and completely DReer vane The 


phase variable representation (Chen, 1970) of such a system 


1S% 
x(k+1) = Ax(k) + bu(k) + w(k) (25Fa)) 
Vak) v=. isk 6+ We Ck} Corer) 
where x(kK)iis an nxile state veotor, 


u(k) is an input variable, 
w(k) is an nx1 system noise, 
y(k) is an output variable, 


y(k) iS a meaSurement noise, 


I 
A = ja|— 

0 
hs= A[ #i0%- @] *¥ 
au= [a, aot; 
bee hbat aera 


I is the unity matrix of appropriate dimension. 


and ' denotes tranpose of a matrix. 


Here w(k) and v(k) are assumed to be N(0,Q) and N(0,r), 


respectively. By means of Kalman filtering theory (Kalman, 
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1960) the number of the unknown parameters which 
characterize noise can be reduced by an innovations 
representation (Kailath, 1968; Geesey, et al., 1969; Astrom, 
1970) to a form which requires a total of only 3n+1 unknown 
parameters to completely characterize the n-th order system 
to be identified against [2n+1]+[n(n+1)/2] unknown 


parameters in equation (2.1), namely 


ie cure e euch: avin Ol (2.2a) 
Hie eae ee (2.2b) 
where ViCkK)O TS NCOUw)S 


w= hiPhi+ r 
P is the steady-state covariance matrix of the 
State error of the Kalman filter, 


and d=[d, -:- d,]' is the steady state Kalman filter 


gain. 


Since the system matrix A in equation (2.2a) contains the 
unknown vector a, it is useful to replace the partially 
known matrix A with F, a known matrix but preserving the 
Same structure of the matrix A for the identification 


problem. 


Equation (2.2) can be thus rewritten as: 
xii eaerx(k) ets Va-t)xqtk)) + buvk) Rady (Ck) (22a) 


y(k) = h'x(k) + v(k) [2.2b] 
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Using equation (2.2b), equation (2.3a) becomes 
mkt?) MSBP x(k) MPCE=6) 0 CR) ee Ca- By Cee e bulk) G2r3:b) 
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y(k) 


where c’=a-d 


Regarding the state vector x(k) as the sum of the 
State vectors x‘')(k), x°??(k) and x‘?)(k) of the three 


Systems driven by v(k), y(k) and u(k), respectively, we 


obtain: 
cus Gkee )) Mem Exe a neem OAKS) | (2.4a) 
x GPECK+ 1) MECPxOIMER) OF (a=pyy(k) (2.4b) 
x63) (k+1) = Fx{?)(k) + bu(k) (2.4c) 
yCkv= Chix Oe) (kK) SHCh eee Ck StS" Kew (Sey CK) (2.4d) 


Since equations (2.4) can be considered as a combination of 
three singleinput singleoutput (SISO) systems they are 


rewritten as: 


mack+ 1), = Fz, (ki, = avck) (25-5a) 
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which are input-output equivalent to equation (2.2). The 
block diagram representation of this system is provided in 
Fig. 2-1 where q denotes the shift operator 


ee eC ey Ky mary: (KH): 5 


This representation has long been used in the adaptive 
observer problem for deterministic systems (Anderson, 1974; 
Kraft, 1976; Kreisselmeier, 1977) as is the case v(k) = 0. 
The apparent advantage of this representation for the 
deterministic system is the linearity in parameters (Astrom, 
et al., 1971) so that the usual least-squares method can be 
used for identification problems (Kudva, et al., 1974). 
However, for stochastic systems this advantage quickly 
disappears since equation (2.5d) is nonlinear in the 
parameters because v(k), therefore z,(k), is not available 
for measurement. It should be noted that the above 
representation is equivalent to an ARMA-model if all the 


elements of f are zero. 


Besides the determination of the structure of model an 
important requirement in the identification problem is the 
choice of an appropriate model order. It is well known that 


the parametric models result in estimates with large errors 
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when the order of models does not match with that of the 
system to be identified. There are several methods 
available for choosing an appropriate model order (Anderson, 
1962: van den Boom, et al., 1974: Unbehauen, et al., 1974; 
Soderstrom, 977+ Young, let al.4. 1980)2) "Since thismsubject 
also stands alone as an interesting research project we 


leave it here without further discussion. 


Since the oe aia of the Sone input u(k) does not 
affect the characteristics of equations (2.5) except to 
decrease the number of unknown parameters, we shall drop the 
Sat input u(k) in the sequel for the sake of simplicity, 


giving rise to a problem-of time-series identification. 


The formation of a state vector 6(k) by including all 


the unmeasurables in equations (2.5) as: 
O@(k)' = (60 ake heels Cis 76 Suek nd 


6,(k) = 2,(k) 
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lends itself to a nonlinear equation in the state space as: 


6(k+1) = ©6(k) + nv(k) (C2.6a)) 
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and O'S are understood as vectors or matrices of 


appropriate dimension. 


2.3 Approaches to the Linear Identification Problem 

Most of the approaches to the recursive identification 
problem of linear stochastic systems have been developed 
from statistical considerations and many of them are based 
on the modification of the LS method. Among them the best 
known are the RML method (Soderstrom, 1973), the recursive 
instrumental variable (IV) method (Young, 1974), the 
recursive prediction error (PE) method (Moore, et al., 1979: 
1980) and the EKF algorithm (Kopp, et al., 1963; Cox, 1964). 
The main shortcoming of these methods, though some progress 
has been recently made (Ljung, 1977b), is the difficulty 
encountered when attempts to analyze the convergence 


properties of the corresponding algorithms are made. 


Recursive identification problems could also be 
treated as a model tracking problem. It then becomes 
natural to introduce the model reference adaptive (MRA) 


method (Landau, 1974) as a design tool. In this approach 
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the recursive algorithm is designed to make the function of 
errorssdecreasesasymptotically: ‘For)thisetaskyhiteissuseful 
to apply well-known stability theories, notably the 
Kalman-Yakubovich (K-Y) lemma (Hitz, et al., 1969) and 
hyperstability theorem (Landau, 1969; 1979). The advantage 
of this approach is that the overall stability of the 
resulting error system with respect to the convergence of 


the algorithm is assured. 


Two configurations related to the MRA technique could 
be laid out depending on the error measurement methods 
(Lion, 1966). The equation error methodt has been widely 
used for the problem of adaptive observers and schemes which 
work satisfactorily in real applications have been reported 
(Luder,Ret-al.» 1974).q@ The poputapigy of this method for 
the particular problem is largely due to the ease with which 
the gradient of the quadratic cost function of the errors 
can be obtained and the simplicity of designing an recursive 
algorithm to decrease the cost function asymptotically. The 
results of this method, however, are known to be biased if 
correlated noise is present. For this reason, the 
applicability of this method is somewhat limited for 
stochastic problems. The output error method, on the other 


hand, has an appealing potential in that it has been shown 


tThe equation error method is often referred as the 
series-parallel approach and the output error method the — 
parallel approach. These terms apparently have their origin 
in the configurations used for error measurement 
(Martin-Sanchez, 1976). 
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to yield an unbiased estimate even in the presence of 
correlated noise (Landau, 1976; 1978). Research has since 
concentrated on the possible extension of this method to the 
problem of the identification (Dugard, et al., 1980) as well 
as to adaptive control of stochastic systems. However, for 
stochastic systems it turns out that the conditions imposed 
for the overall Soa ey are not always possible to be 
Satisfied because of the involvement of the parameters that 
are to be identified (Ljung, 1977a; Egardt, 1980; Johnson, 


iS80). 


Fig. 2-2 shows the block diagrams for two error 


measurement methods 
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The inclusion of the unknown parameters to the 
augmented state vector transforms a parameter identification 
problem into a nonlinear eStimation problem. Since the 
optimal solution to nonbinear-filterung problems 1s not, in 
general, realizable in the finite dimensional systems (Bucy, 
et al., 1971), nonlinear filters-are, in most cases, 
obtained by approximate design methods. While the EKF 
offers the best simplicity in implementation, thus, being 
most widely used among the numerous methods, the application 
of this approach often unnecessarily increases the 
computational requirements if only the parameters need to be 


known. It would appear to be worth while to attempt a 
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a) The equation error measurement 


b) The output error measurement 


2.2 Two error measurement methods 
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modification of the EKF in order to minimize the 
‘computational burden. It is noteworthy that the convergence 
properties of some modified algorithms of the EKF could be 
analyzed with the stability theories used in the case of the 


MRA method. 


In the néxt chapter, we will first discuss the 
identification problem for the deterministic system and 
develop a fast algorithm through a clear insightful 


understanding of the existing algorithm. 
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CHAPTER: III 
A RAPID IDENTIFICATION SCHEME FOR DETERMINISTIC SYSTEMS 


3a Vel ntroducticn 

Since the stochastic identification problem formulated 
in the previous chapter is intrinsically nonlinear and thus 
Practically not amenable to an exact solution, an 
approximate solution will be sought through the 
linearization of the problem. While the setting in which 
the solutions sought for the deterministic identification 
problem can be considered to be too idealistic from the. 
stochastic point of view, the technique developed in this 
chapter for the improvement of an existing deterministic 
identification algorithm will be useful to conceptually 
understand and thus solve the stochastic problem. The 
remaining portion of this chapter will be devoted to the 
application of MRA techniques to the deterministic 


identification problem. 


As briefly mentioned in the previous chapter, the 
preservation of the linearity in the parameters, by which it 
is meant that a generalized error is linear in parameters 
(Astrom, et al., 1971), of the deterministic system provides 
a good opportunity for a wide range of applications of the 
MRA technique to various problems of the linear 
deterministic system. Furthermore, the linear system 


theories represent the most active research area and thus 
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are readily available for various applications including MRA 
schemes. Consequently, numerous successful MRA schemes for 
the identification problem as well as the control problem 
have appeared in the literature (Landau, 1974; 


Kreisselmeier, 1980). 


With a few exceptions, most of the recent developments 
in the MRA schemes have adopted the stability-based approach 
in designing an adaptation mechanism. -. In this approach, 
Stability theories are used to determine the sufficient 
conditions under which the adaptation mechanism makes the 
overall system asymptotically and globally stable. Thus, 
the design problem associated with the deterministic MRA 
method corresponds to a stability problem. The crucial 
point in the design procedure is therefore to choose an 
equation for updating the adjustable parameters in such a 
way that the overall nonlinear error system, which 
constitues the state error system and the adaptation 


mechanism, is asymptotically stable. 


Kudva and Narendra (1974) employed the Lyapunov Direct 
Method, which is one of the most widely used stability-based 
design methods, to develop a simple adaptation mechanism for 
the identification of the discrete multivariable system. In 
the design procedure, it is assumed that all the state 
variables are accessible and that there is no noise 
involved. While this identification scheme (K-N scheme for 


short) possesses the overall stability, simulation studies 
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show that the convergence speed is often undesirably slow. 


Thus, there is need for increasing the speed of convergence. 


A geometric interpretation of this adaptive scheme 
provides a greater insight into and consequently a clearer 
understanding of the scheme. Though the underlying concept 
1s rather simple, this becomes very crucial in the | 
development of an algorithm with rapid convergence property 
which is important for on-line applications. This approach 
also provides a general idea of how the gain of the improved 
algorithm is to be modified if noise is present. This, in 
turn, has an enormous implication in the development of an 
identification algorithm for the stochastic system developed 


in Chapter IV. 


3.2 Geometric Interpretation of the K-N Scheme 

Kudva and Narendra (1974) considered the problem of 
identifying the unknown but constant nxn and nxr, (n2r) 
matrices A and B of a noise free dynamic discrete system 


given by 
x(kt+1) = Ax(k) + Buk) Coty) 
and have derived the following adaptive algorithm: 


Peke1) = C8(k) + [AC k+1) - C]x(k) + B(k+1)u(k) (3.3) 
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PRC ome (cates (ACh) BCR 


axfe(k) - Ce(k-1)]@(k-1)' 
x fe Se eee ee eames C373) 
AO(k-1)'d(k-1) 


where C is an nxn stable matrix, i.e. all its eigenvalues 


lie inside the unit circle, 

Oe key Meares) Meeks) Ie 

A 1S the largest eigenvalue of Z 
and rea CkYy ex Ck) 


Kudva and Narendra have shown that estimate 
[A(k),B(k)] will converge asymptotically to [A,B] provided 
that 

i) the gain matrix Z is chosen to be symmetric and 

positive definite, 
Hie)A0 a< Fa KS 25 
iii) the system (3.1) is completely controllable 
and iv) the input sequence u(k) is general enough (Kudva, 
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where P(k) = B(k) - P, 


P(k) = [A(k),B(k)] 


and Pr=rA, Bs. 


Since the matrix C does not appear in equation (3.5), the 
convergence of this equation which is the main concern here, 
does not depend on whether C = 0 or # 0. In simulation 


Studies, C = 0 is the most convenient choice. 


In proving the stability of equation (3.5), for the 
Sake of convenience, we institute the following change of 


variables: 
w(k) = @(k-1)/] | o(k-1) | | (3.6) 


where ||.|| denotes the Euclidean norm. 


Equation (3.5) reduces to: 
P(k+1) = B(k) - aZP(k)W(k)w(k)'/d (35.475) 


Since. P is a positive, definite, matrix, there exists a 


nonsingular matrix T such that 
ToT GSN (3:28) 


where A is a diagonal matrix whose diagonal elements A; lie 
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Premultiplying equation (3.47)) bys? vand sing equation 


(3.8), we get: 


P(k+1) = P(k) - aAP(k)W(k)W(k)! (3.0) 


where P(k) = TB(k) eae 


We can partition equation (3.10) row-wise as below: 


Bi (k+1) = 6 (k) - ad W(k)W(k)'D; (k), ieft,n] (Gna) 


where p,(k)' is the i-th row vector of the nx(n+r) matrix 


P(k), A; is the i-th diagonal element of A 


and OM Tae eGo 
Since 

pi(k)'wtk) = ||pi(k)||cosé,, (ge) 
the choice ad; = 1 makes ||6;(k+1)|| minimal for a given 


w(k) where 6, is the angle between p,(k) and w(k). 


If we choose a = 1 for convenience, A, must be 1 and 


A = 1. Furthermore, if we choose T = I for convenience 
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3.1 The geometric effect of gain 


Fig. 3-1 shows the geometric gamennpareciion of 
equation (3.11). Now that we have obtained the optimal gain 
aZ/d = I, we can go back to the original equation (3.7). 
Wiene ties chovce sof 2 = land a°=, 1, equation (6.7) cans be 


written as: 
Di (k+1) = B(k) - W(k)w(k)'B, (k), iel1,n] (m3) 


where p;(k)' is the i-th row vector of P(k). 


We can see that for a given wW(k) the best possible scalar 


Gain atv each-step isiW(k)palk) > 


The global asymptotic stability of equation (3.13) is 
guaranteed by the fact that condition (iv) in the K-N scheme 
implies that any nt+r consecutive vectors ¢(k) span the ntr 


dimensional Euclidean space, i.e. 
botk). O(k#1)) 4. eek ean toa) is nonsingular. 


The implication of equation (3.13) is shown geometrically in 


Figure 3-2 for third order $(k)'s. 


Keeping in mind ‘the directions of $(k-1) and W(k) are 
the same, observation of Fiure 3-2 reveals the following 


properties of pj (k): 


ave p, (ke 1)) 1snenenogonal™ to W(k), 


b) p,(k+1) lies in the plane made up of p,(k) and w(k) 
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c) if pi(k) is parallel to w(k) then 5,(k+1) is zero, 


and d) ||/Bi(k+1){| < | 1B; (k)| 


e 


3,3 A Rapidly Convergent Algorithm (Kim, et al., 1987) 

While equation (3.13) maximizes the decrease of the 
parameter error p,(k) at each step and thus p;(k) is assured 
of convergence to the origin, the arbitrariness in the 
direction of w(k) retards the speed of convergence 6,(k) to 
PhewOrigin. in Order to circumvent this dit ficubiys awe 
introduce an orthogonal sequence of ¢(k) in the place of 
Ww(k) to coordinate systematically the direction in which 
pi (k) decreases and adjust the scalar gain a(k) in sucha 
way as to make p;(k+1) orthogonal to W(k) (This adjustable 
gain a(k) should not be confused with the fixed parameter a 
used in Bhai 3.2). srhen it 1S!nGhidviticultero see from 
the properties a) and b) that this n+r vector can be made to 
be parallel to v(k) in ntr-) steps. “lpothis is done, then 
by property c), p;(k) will converge to the origin in one 
more Step... Thus the  totali number of steps required to bring 


Bi(k) to the origin will be utmost ntr. 


Another way of explaining this is as follows. We 
choose w(0) as the first basis vector of the nt+r Euclidean 
space, v(0) = w(0). We select W(k) as the remaining basis 
vectors. We also choose a(k) to maximize 
[1b:(k)|| - | ]pi(k+1) || at each step so that the number of 
steps required to bring 5;(k) to the origin will be a 


minimum. 
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In the light of the above discussion, we can use the 


following algorithm in place of equation (3.13). 
Bi(k+1) = Bi (k) - alk)W(k)o(k-1)'B,(k), deL1,n] 9 (3.14) 
The) choicesofia(k) to maximize | 15; (k)|1*-9]|8; (k+1)]]| is 


made as follows. 


If an orthogonal sequence W(k) is generated from the 


sequence of ¢(k) by the Gram-Schmidt process, the plane made 


up of ¥(k) and 5;(k), and the plane made up of W(k) and 


@(k-1) are orthogonal to each other (see Figure 3-3). p,(k) 


can be expressed as: 


Bi(k) = BPiylk) + Bry(k) 
where B,y(k) is a projection of 6;,(k) on W(k) and Daa Ck us 
a projection of 6,;(k) on a vector orthogonal to both v(k) 
and @(k-1). Since @(k-1) is orthogonal to Biy(k) and 


v(k)'W(k) = 1 because of orthonormality, we can write 


Olk—t) pith) sae k= puke) 


g(k-1) "(0 (k) 0! (kK) DB; (kK) ] 


[o(k-1)'O(k) JID (Kk) 'B; (k)] (2s) 


Using equation (3.15), we can write equation (3.14) as: 
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Bi (k+1) = Bi (k) - alk) w(k)(o(k-1) (kK) ]0(k)'B; (k) (3. 16) 


fore eli -0 | 


Following equation (3.13), it has been pointed out that 
w(k)'B,(k) is the best possible scalar gain for a given 


w(k). By analogy, we see that the choice of 
a(k) = 1/[W(k)'@(k-1)] (3.17) 
in equation (3.16) will make (k)'D,(k) the best possible 


gain for a given W(k). 


The resulting equation: 


0(k)¢(k-1)'B; (k) 
pi(k+1) = 6, (k) - —e—__———— _iel1,n] (3.168) 
o(k-1)'W(k) 


is basically the same as the equation (3.13) except that the 
direction of the sequence of 6, (k+1) - p;(k) is now 
orthogonal to that of all preceeding elements, which enables 


pi(k) to converge to the origin in utmost n+r steps. 


The composite equation for all i, ie[1,n] would be: 


P(k)o(k-1)W(k)' 
(3.19) 


P(k+1) = P(k) - = 
ye CK— 1) 
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Thus the resulting algorithm is: 


e(k)(k)! 
Bake erp (ky ope (B20) 
W(k)'o(k-1) 


Figure 3-3 shows the goemetric implication of equation 


Va-ake), 


The algorithm proposed by Westphal (1978) yields the 
same results as the, algorithm (3.20)wonw ft) pasinaciie 
Westphal's scheme is set equal to n+r-1. However, the 
emphasis placed by Westphal on the utilization of a priori 
information about the system parameters to improve the the 
convergence unnecessarily obscures the real strength of the 
algorithm. As the simulation results show in the next 
section, a priori information is not at all necessary for 


ensuring rapid convergence of the algorithm. 


This algorithm is also readily applicable to a 
deterministic MIMO system whose present output can be 
described as an inner product of the parameters and a 


sequence of the past outputs and inputs. 


3.4 Simulations 


Using the algorithm (3.20), a computer simulation is 
first carried out to identify the parameters of the 


dynamical discrete system: 


x(k+1) = Ax(k) + Bu(k) 
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3.3 The geometric effect of orthonormal vectors 


Te Nie 


fi 


erosaev Ismionoitao Fee debts a 
" iT O : : : 


: | | iw oo | 


where A and B are the unknown but constant matrices. For 
the purpose of simulation, a fourth order discrete system 


with the following values for A and B matrices: 


Pe Oe Baal 2. 6990R=4 I; 64575-4494) 555425—3 
4,7943E-4 SO IIR) Wa GOOb a a OC Or ag 
F991 IEes 3.6498E-3 FS USE! 1.4074E-2 


5.0006E-6 hacsU Tea 5 92 JOOUnas UO On oes0 


a a A lo i. /O43E=3 
Soe 2 ha Zeeman OOS Zio 
-5.4940E-2 4.4605E-2 


Zee OA 2.230 \E=8 


is chosen to demonstrate the effectiveness of the algorithm. 
This system is actually obtained by discretizing the 
continuous system in (Narendra, et al., 1973; Luders, et al, 
1974) with sampling period 0.01s. Here the initial values 
of A(k) and B(k) are chosen to be all zero and the matrix C 
in equation (3.2) is conveniently set to be zero. A part of 


the results of the simulation is shown in Figure 3-4, 


Another simulation is carried out in order to show the 
usefulness of the proposed algorithm to tackle a more 
realistic problem which was generated by assumption that 


only one state variable is measurable in the above example. 


mck { a= OAs UK) + Buy) 
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y(k) = [0 0 0 10]x(k) 


This system is then equivalent to the system in equation 


(2.1) in terms of the input-output relation as: 


SCkt 1) eae Ax Ck) +ibpuy(k) + bous Ck) 


y(k) = h'x(k) 


where 

I 
A = ja|—— 

0 
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Following the procedure discussed in Section 2.2 to 
represent the above system in a more convenient form for the 


identification problem, we obtain: 
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y (k) = ez ke a byezeek) 22 beteky) 
where F, is described in equation (2.3) with f = 0. 


Then, we can readily apply the identification algorithm 
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proposed in Section 3.3. 
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CHAPTER IV 
ALGORITHMS FOR IDENTIFICATION OF STOCHASTIC SYSTEMS 


4.1 Introduction 

In this chapter we shall develop some new algorithms 
for the identification of the stochastic system described in 
ChapterseI,d. jiSince: thes system! dsestochastic) 4 tedisealso 
natural that we consider the use of stochastic estimation 
techniques in this development. It was pointed out in that 
chapter that the augmentation of the state vector of the 
stochastic system with the unknown parameters makes the 
error system nonlinear and consequently we have a nonlinear 


estimation problem of a stochastic system. 


The main objective of stochastic estimation problems 
is to obtain the conditional mean values and covariance 
matrix of the unknown quantities. The mean is the estimate 
with the minimum variance for the given measured data and 
the covariance matrix represents the measure of the 
uncertainty in the estimates. It is well known that the 
conditional mean and covariance matrix cannot be, in 


general, characterized by a finite set of the moments of the 


conditional density function. 


The involvement of solutions of functional integral 
difference equations (Bucy, et al., 1968) makes the 
computation of the exact conditional density function, in 


general, virtually impossible. Thus, in practice, 
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approximations of the conditional density function are used 
and this results in suboptimal filters. it is rather 
desirable that the conditional density function during the 
process of approximation is represented by a finite set of 
Parameters so that the corresponding nonlinear filter is 
made up of finite number of equations of evolution for these 
parameters. Unfortunately this simplification cannot be 
achieved without the loss of sufficient statistics on the 


estimates. 


It is known that in the linear estimation problem the 
conditional density function, being Gaussian, can completely 
determine the mean vector and covariance matrix. In this 
case, the mean vector is uSually treated as the state of a 
linear filter. This filter, in turn, can be looked upon as 
a reference model in a model tracking context. The 
covariance matrix is used to compute the direction and 
magnitude of the gain vector. In order to take advantage of 
this useful feature of the linear estimation problem, the 
nonlinear estimation problem is linearized by taking only 
the hiinear spart i6fiithe tTay.lor sertesipepresentaaion io& the 
nonlinear system. If this ys done, then, the Kalman 
estimation theory can be readily applied to the approximated 
system because at each step only linear terms are involved. 
The resulting filter is called the Extended Kalman filter 
(EKF), which has been most widely used in real-life 


applications (Jazwinski, 1970). 
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It 1s known that the convergence of the estimates of 
the EKkPealgorithm to, the. true) values® istassured iftithe 
estimates stay within a small neighborhood of the true 
values. In the case of the MISO system described in Chapter 
II, the augmented state vector consists of the state vector 
and the parameter vector. If the parameter estimates stay 
within a small neighborhood of the true values, generation 
of good estimates of the state vector automatically follows 
in the innovations representation. There is no need to 
estimate the state vector separately. However, use of the 
Standard EKF algorithm requires that the estimates of the 
entire augmented state vector be determined at each step. 
This unnecessarily increases the computational burden. One 
of objectives in this chapter is to develop certain 
modifications of the EKF algorithm to reduce this 


computational burden to facilitate on-line applications. 


Stability is perhaps the most important consideration 
in assesing the performance of any techniques of recursive 
adaptive schemes. In the deterministic system, a number of 
methods are available to test or analyze the stability of 
recursive adaptive techniques, for example, Lyapunov method 
and the hyperstability criterion. Among these, the 
hyperstability criterion, which can be especially regarded 
as a specific formulation of the Lyapunov method, offers the 
most systematic approach to system design. The application 


of the hyperstability criterion to a stochastic 
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identification problem was first pioneered by Landau (1976) 
and became a popular tool to design a recursive adaptive 


scheme. 


In Section 4!..2 the application of the EKF algorithm to 
the identification problem is discussed. In the remaining 
eections Of the chapter, several modified /eKE algorithms for 
the stochastic identification problem are developed and 
certain similarities of these algorithms with earlier 
algorithms (Panuska, 1968: 1969; Soderstrom, 1973) are 


shown. 


4.2 EKF Algorithm for Identification Problem 

In this section, our aim is to derive EKF equations 
for the identification of the system described in Chapter 
II. The derivation follows the approach in (Anderson, 


Stealer oO 79:)i, 
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and O's are understood as vectors or matrices of an 


appropriate dimension. 


The nonlinear function, which appears only in the 
observation system (2.6b), can be expanded in Taylor series 
as: 


y(k) = y(k) + W(k)'O(k) + 1/26(k)'SO(k) + v(k) (Araie} 


where $(k) = 6,(k)'6,(k) + 22(k)'63(k) Cae.) 
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of 6(k) based on the observations up to k-1 step. 


Neglecting the second order term in equation (4.1), we 


obtain a linearized approximate version of the system (2.6) 


as: 
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O(k+1) = 60(k) + nv(k) (4.5a) 


y(k) = $(k)'O(k) + v(k) + $(k) 7 (4.5b) 


Then, the application of the Kalman estimation method to 


equation (4.5) yields the following algorithm: | 


B(k+1) = 86(k) + glk)Ly(k) - §(K)] (4.6a) 
y(k) = 6,(k)"62(k) + 22(k)'63(k) (4.6b) 
g(k) = [@P(k)W(k) + nw]m(k) (4.6c) 
P(k+1) = @P(k)S' - g(k)m(k)~*g(k)' + nwn' (4.6a) 


P(0) = -P(0)" 2-0 
m(k) = [w(k)'P(k)w(k) + w)-' (4.6e) 
where w is the covariance of v(k) and P(k) is the 
approximate covariance matrix of EKF. 
g(k) and P(k) can be also, after a matrix manipulation (see 
APPENDIX 4A), expressed as: 
g(k) = 6,(k)P(k)w(k)m(k) + 7 au) 


P(k+1) = @,(k)(P(k) - P(k)w(k)m(k)W(k)'P(k)]&,(k)' (4.8) 


where ©,(k) = ® - nw(k)' (4.9) 
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case g(k) approaches 7. 


The error state @6@(k) satisfies the difference 


equation: 


@(k+1) = 6(k) - g(k)Ly(k) - ¥(k)] + nv(k) 
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- ,(k)P(k)W(k)m(k)v(k) (4.10) 


As clearly seen in equation (4.10), the effect of the 
nonlinear term 6(k)'S@(k) is obvious. The convergence of 
the error state 6(k) to the origin could be thus assured, 
only in a neighborhood af he origin of the error state 


Space where the effect of thenonlinearity is negligible. 


Another advantage of employing the innovations 
representation, besides the reduction in number of 


parameters to characterize noise, is that the solution to a 
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in the neighborhood of the origin of the parameter error 


Space, where the EKF is assured of convergence to a true 


values of the parameters. 


4.3 First Modification of the EKF Algorithm (EKF-M1) 
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Suppose we assume that 6,(k) is nearly equal to 6,(k), 
then P,,(k) which is an approximate measurement of the 
uncertainity in 6,(k), can be set equal to 0, thus reducing 
the computational requirements. A convenient choice f = 0 
considerably simplifies P;2(k) and P22(k) for computational 
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It should be mentioned that a certain numerical problem 
associated with equations (4.16-17) will arise such that 
P22(k), which approximately represents the uncertainity in 
the estimates 6.(k) and 63;(k), could be nonpositive 
definite. In order to get around this difficulty, we simply 
add a positive scalar quantity uw to m,(k)~' so that P22(k) 
remains positive definite for all k. Then this algorithm 
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m,(k)(U(k)'6.(k) + $(k)]'Po2(k) (A202) 


Since in a small neighborhood of the true values of the 


parameter in [62',@3;']' vector space, it can be assumed that 
Poo(kt+1)P22.(k)-! & 1 (4,23) 


Hence we can write U(k) and P22(k) as: 


U(k+1) = [Fo' - hO,(k)'JU(k) - ho(k)' (424) 
Peek yin = eR os (kh) =P o bk OCH i nm. Uk) hotCk area ek) 
(4.25) 


If we choose a positive scalar quantity u such that 

w= 62(k)'U(k)P22(k)U(k)' 62 (k) (4.26) 
then we obtain 

m,(k) = [h'U(k+1)P22(k)U(k+1)'h + w)! Gen2%) 


so that P22(k) is assured to be positive definite for all k. 


Hence the modified EKF becomes: 
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Po2(k) and U(k) are described in equations (4.24-25). 


The salient feature of this algorithm, which we call the 
first modified EKF (EKF-M1) in the sequel, is the retention 
of P,2,(k) so as to preserve the local convergence property 
which is considered to be a very useful feature of EKF. 
However, Since this modification of EKF is based on the 
assumption that the nonlinearity effect is negligible ina 
neighborhood of the origin of the parameter error space, it, 
in general, suffers from the lack of the global convergence 
property. However, under certain conditions, such global 
convergence is possible. The convergence of P(k), as shown 
in equation (4.8), completely depends on the stability of 
@(k)) which, in turn; is determined by 6,(k). Lesis 
therefore necessary to ensure that §.(k) lies in the region 
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restrict the extent of observation nonlinearity. A recent 
study (Ljung, 1979) has shown that the modified EKF 
algorithm (4.28) when complemented by a facility to confine 
6.(k) to the region where F,' - hé,(k)' is stable, will 
enable the parameter estimates to converge to the true 


values with probability 1 (w.p.1). 


It will now be shown that the approximate RML 
algorithm (Soderstrom, 1973) can be Ber es from EKF-M1 as 
follows. Since 6,(k) is approximately equal to 6,(k) ina 
neighborhood of the parameter error space, it follows that 
Sukie iy Gees iy Ck) cis equal.to the wnoise term vik)? This: 
in turn, suggests that 6, (k) can be obtained from equation 
G2,.5a) using e(k)) instead of v(k).. If.this 1s done; the 
algorithm turns .out to be the same as the approximate RML 
algorithm with w = 1. That the approximate RML method 
possesses the local convergence property (Ljung, 1975) can 
be seen by noting that the gain gz(k) in equation (4.28e) is 
almost equal to h in a neighborhood of the origin of the 
parameter error space. However, the global convergence 


remains in question (Holst, 1977). 


4.4 Second Modification of the EKF Algorithm (EKF-M2) 


P,.(k) represents the correlation between two error 


vectors, 6,(k) and p(k) where 
6.(k) = 0,(k) - 64(k) (4.29) 


B(k) = p(k) - plk) (4.30) 
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The assumption that 6,(k) = 0, however, reduces the 
importance of P,2(k), thus providing grounds for further 
Simplification by setting P,2(k) = 0 (thus U(k)) in equation 
Care2 8) 


This results in the following algorithm: 


SET VND oncns etre striae tar iG 33228) 
6,(k+1) = Fo'6,(k) + hmz(k)ly(k) - $(k)] (4.32b) 
y(k) = o(k)'p(k) (4.32c) 
g3(k) = Po2(k)o(k)m2(k) (4.324) 
Poo(kt+1) = Po2(k) - Po2(k)o(k)m2(k)@(k)'P22(k) (4,.32e) 
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In the sequel we shall refer to this algorithm as the second 


modified EKF algorithm (EKF-M2). 


Following the same reasoning as before namely e(k) 
could be assumed to be equal to v(k), we use equation (2.5a) 
with e(k) in the place of v(k) to obtain 6,(k). Then the 
resulting algorithm with w = 1 in equation (4.32f) turns out 
to be the ELS method derived by Panuska (1968; 1969). The 


main drawback of the ELS method as well as the second 
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modified EKF algorithm lies in the requirement of additional 
restrictions on the system in order to ensure the 


convergence of the estimates to the true values. 


4.5 The Convergence Property of the EKF-M2 Algorithm 

In this section, the convergence property of the 
second modified EKF algorithm will be discussed using the 
Kalman-Yakubovich theorem which up to now has been used to 
study the stability of deterministic problems. Recall the 
brief discussion in section 2.3 that the identification 
problem can be formulated as a model tracking problem. In 
order to analyze the convergence property of EKF-M2 in the 
model tracking context, we shall examine the stability of 


the error model. We obtain the error states 6,(k) and p(k) 


satisfying the following equations with assumption w = 1 as: 
6,(k+1) = Fo'6,(k) - hm3(k)e(k) + hv(k) (4.33a) 
pok+t). = p(kie=4g Ck) Gk) (4.33b) 
e(k) = y(k) - ¥(k) 
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obtain the following feedback system: 


6,(k+1) = [Fo’ - h62']6,(k) - hy(k) (4.36a) 
ex(k) = 62°6,(k) + y(k)+ v(k) (4.36b) 
Dik ye= sp Uk)ij7 Peek) (kK) eq. Ck) (4.36c) 
y(k) = o(k)'p(k) - $(k)'P22(k)o(k)e,(k) (4.36d) 
Poo(k+1) = Po2(k) - Poz2(k)o(k)m3(k)¢(k)'P22(k) (4.36e) 


P22(0) = P.-(0)' AG 
Since y(k) can be expressed as $(k)'D(k+1), we get 
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Comparing the above equation with (4.33c) we can easily 
notice that e,(k) is the error after a correction is made on 
B(k). In the sequel e,(k) is referred to as the a 


posteriori error while e(k) is the a priori error. The 
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condition for the convergence of the EKF-M2 algorithm to the 


true values will now be stated as below. 


theorem: is Suppeses that Fou- hé>" tin equation (436a) has 
all its eigenvalues inside the unit circle, then the errors 


6,(k) and B(k) satisfying equation (4.36) vanish as k > @ 
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is strictly positive real (s.p.r.) 


Wheres crt = fey ice] 


and q is the shift operator. 


The proof of Theorem 1 consists of two steps. First, we 
will prove using the Kalman-Yakubovich theorem that 6,(k) 
and p(k) converge to the origin without v(k) in equation 
(4.36b). Then we prove that v(k) actually has no effect on 


the convergence of 6,(k) and p(k) to the origin. 
We need the following lemma. 


Lemma: Given the following set of time-varying systems: 


Be tk+d) = FO, Ck) = hy Gk) (4.39a) 
e,(k) = 02'8,(k) + 7(k) (4.39b) 


B(k+1) = Blk) - P22(k)o(k)e,(k) (4.39c) 
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Oe OCR) Dba 6k) "Pes (kK) o Ch emck) (4.39d) 
Po2(k+1) = Po2(k) - Po2(k)o(k)m3(k)O(k)'P22(k) (4.39e) 


P,2(0) = P22(0)' et) 
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the errors 6,(k) and p(k) converge to the origin as k > @ 


asf BV Gase) Pam oe ees OT iia emer (4,40) 


where F, = Fo - @6zh' 1S assumed stable 


proof: See APPENDIX 4B. 


The above lemma is the discrete version of the 
Kudva-Narendra theorem (Lin, et al., 1978) with a 
modification that the gain matrix P22(k) decreases in 
magnitude as k increases. The decrease of P22(k) is 
necessary to account for the rejection of the norse in 
equation (4.36b) by gradually decreasing the weight on the 
noisy observations. It, however, requires a more stringent 
restriction on the system than is neccesary with a fixed 


PoAGk) GLandau, 1979). 


Resuming the proof of Theorem 1, note that the real 
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F matrix. 


From equations (4.34-35) we have: 


ex(k) = [62'6,(k) + o(k)'B(k) + v(k) Jm3(k) 


(4.41) 


v(k) = {@(k)"B(k) - 6(k)'P22(k)O(k) [62°68 ,(k) + v(k) I}. 


ms, (k) 


Substituting (4.41-42) into (4.36), we have: 


KOREN = Ck) xk) + BC k)v(k) 


where 
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Hence, x(k) can be expressed as: 


(4.42) 


(4.43) 


(4.45a) 


(4.45b) 


(4.45c) 


(4.45d) 


(4,46) 


(4.47) 
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col aa K-2 k-1 

x(k) = ME(i)x(0) + EM E(i)8(j)¥(3) 
= 1=7 


=O 5 =j+1 


teks i vCka) | (4.48) 


Since =(k) is stable from Lemma and v(k) is a white noise, 


ie a Giese 0 (4,49) 
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where E(-) denotes expectation. 


Thus (Mendel, 1973), 
Lime Bix(k)} = 0 (47,52) 
kro : 


we have proved Theorem 1 which, in turn, establishes the 


convergence of the EKF-M2 algorithm to the true values. 


The difference between the ELS algorithm and EKF-M2 
algorithm (4.32) is that in ELS a priori error e(k) is used 
for the generation of 6,(k) while the EKF-M2 algorithm 
(4.32) uses a postériori error e,(k). Simulation studies 


show that this enables the EKF-M2 algorithm to converge more 
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rapidly than the ELS algorithm. 


A limitation of the EKF-M2 algorithm, however, lies in 
the fact that the positive realness requirement (4.40) for 
unbiased convergence depends on the unknown parameter vector 
62. A method which helps to overcome this limitation to 


some extent 1s discussed in the next section. 


4.6 Third modification of the EKF (EKF-M3) 

In the modifications of the EKF algorithm discussed so 
far, we have set the vector f in the matrix F equal to 0. 
In this Section, tHe ‘implication jof using a general f {mot 
neccessarily 0, will be discussed, The vector f£ will be 
chosen such that a stable polynominal of the difference 
between f and certain system parameter satisfies a positive 


realness condition. 


If a general vector f is introduced in equation 


(4.32), the following algorithm is obtained: 


Z2(k+1) Fz (kK) aby Ck) (4,.53a) 
O (kt?) =F 16, (k) tehme (kby (kh) 29 Ck) (4.53b) 
B(k+1) = P(k) + P22(k)o(k)ms(k)Ly(k) - ¥(k)] (4.53c) 
y(k) = o(k)'B(k) (4,53) 
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This algorithm is referred to in the sequel as EKF-M3. 


Convergence of the EKF-M3 algorithm 


Recall from equation (2.5) that the system can be described 


as: 
Ora ktas) = bee wCke): +a hy (Ck) ; (4.54a) 
Zeke yee Fs Zak) re Hiyotky) (4.54b) 
Vikas WO Ok) woo ky pet vik) (4.54c) 
where 12) = [0:8 OST" ae (4.55) 
62 = fe 
63 = a-f 


We obtain the following set of error equations: 


6.(k+1)>= F'6,(k) = hmg(k)e(k) + hv(k) (4.56a) 
B(k+1) = p(k) - Paz(k)o(k)ms3(k)e(k) (4.56b) 
e(k) = y(k) - ¥(k) 
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makes it clear that (4.56) constitues a block feedback 


system-described in the following equations as; 


6,(k+1) = [F' - h62']6,(k) - hy(k) (4,59a) 
ei,(k) = 62'8,(k) + y(k) + v(k) a (4.59b) 
ekg jie (i) - Po2(k)o(k)e,(k) (4,59c) 
yk) = 6(k)"Blk) - 9(k)'Paatk)olk)es(k) (4.594) 
Po2(k+1) = P22(k) - Po2(k)o(k)m,(k)¢(k)'P22(k) (4.59e) 


Pi> (0) = P.2(0)' 2S 0 


Fig. 4-1 shows the block diagram of the above feedback 


system. 
A theorem similar to Theorem 1 can be stated as follows 


Theorem 2: Suppose that F' - h@2z' in equation (4.59a) has 
alll its eigenvalues inside the unit circle, then the errors 


6,(k) and B(k) satisfying equation (4.59) vanish as k > @ 
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proof: The proof is entirely analoguous to that of 


Theorem 1. 


In recent years, the use of the approach utilizing the 
Supermartingale theory to the convergence proof of recursive 
identification algorithms has appeared in the literature 
CSolo, 1979- Goodwin, et. al, 1979, Landau 1987:. 1982)..> In 
this approach, it is necessary to find what is known as a 
corresponding supermartingale for the error feedback system. 
Since a stochastic Lyapunov function in the error feedback 
System should be positive and decrease in magnitude with 
time, it can be looked upon as a supermartingale. It is 
Shown below that the use of this idea in the convergence 
proof for the EKF-M3 algorithm leads to a stronger result, 
namely, martingale convergence compared to the mean-square 


convergence obtained in the earlier proof. 
Proof of Martingale Convergence 


This proof should be read in conjunction with 


APPENDIX 4D where the definition and lemmas are given. 


Since, as stated previously, 6,(k) becomes zero as 
p(k) converges to p as a consequence of the use of 
innovations representation, we only need to prove that p(k) 


actually converges to p if the positive realness condition 


(4.60) holds. 
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Now, let s(k) be defined as: 


s(k) & puke P22 (k)-*D(k) C4563) 


where p(k) = p(k) - p as defined in (4.30) 
Then 
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Taking conditional expectations in (4.62), we obtain 


Eis (Uke) + Zee ye Ck Ls CK) es es CK eo Ck (4. 
wmhererec(k)? = Siy (ka ACk-a) 4 (4. 
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Strictly positive real: system, the inequality 
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always holds (Landau, 1969). 


Defining a new martingale r(k) 
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Bir(k+1){7(k-1)} s r(k) — c(k)/k + o(k)?/k (4.67) 
Since y(k) is the signal of the stable feedback system, 
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It follows from Lemma 5 in APPENDIX 4D that r(k) converges 
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Since y(k) is the output of the assumed stable system and 
e(k) is the signal of the stable feedback system, y(k) and 
e(k) are bounded, so is $(k). Hence T(k) is upper bounded. 


This fact and s(k)/k * 0 lead to the conclusion that 


68 


ah oreo By 


7 


c . 


ir. 


en ee. 7 


ioe A Aiea — ay 


faa.) ) 


- j 


ielae a \ Aha + sich eee a 


nete7? 2sedhes? + icain ne fdngie sis al Oe 
ae a» De 3 
| eae. eee: ‘a. Sh a 
@egzavnon (A)y Dade Ge AIRUIGGA “ot Beowmad mar 
bas (|.a.96) | eee ‘bdadose Aid 4 eidginey mobnez 


(G8 .>) t .giak a> Nays ke 

; Segue +7 ene 

Gre tC sbee avis bang ste Giri | Bas (Ne aot 
.) ee ope ae tine bivode CEI 


ni sian ov (area « tor paitzea 
(ayeda aiaperne 7 
a.e) Calan cremene ~ {AD = (74st 
eee ae ~ SH vie ah og 


ta ja lx wi ate ays 


(— (aS # AS 


bas mecta¢e eidete beupeee @f2 Jo tuqdue sfd ai (Ady sonte 
bas (d4)y ,metdve ssedbee? @idste eft to fenpia sd’ ei Ga)» : 
,babadgod Tagg ei (#17 eoneh .(x)*® 2] of ,Bebaved sus a. ; 
dans noleu Portog oq3. of Ges! 6 + iN dhe Bas 3263 bint 7 


4 Sid aoa all 


Eimte OCK) =O! wpe 1 


k >a 


Since (4.60) is equivalent to (Siljak, 1968; 1970) 
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Pts Corea eee wen a eos (Cc. Eig ia 2 
if the polynomials 
Bc) = |- fiacad re pare by Carer 
Gye.) =a Cua: S665. Cinclane (4927.25) 


are stable, the real strength of Theorem 2 is that the 
estimates 6,(k) and p(k) are assured of convergence to the 
true values kiff'a “vector fi tisiichosen such that -i-c1..e.) 4.64%, 
Satisfies the condition (4.70) regardless of the vector c 
unlike in the case of the EKF-M2 algorithm (4.32). 

Theorem 2 defines the extent of the sensitivity of the 
erconithm. (4.53) to,themae prion knowledge of the 
parameters. In fact, the algorithm (4.32) could be 


considered as a special case of the algorithm (4.53) with 


f = 0. 
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APPENDIX 4A 


[@P(k)W(k) + nw)lW(k)'P(k)W(k) + w]e? 


{(@ - nv(k)'I]P(kK)W(k) + nlwlk)'P(k)W(k) + w]} 
[W(k)'P(k)W(k) + w)-! ; 


[C@ - nw(k)']P(k)W(k)[W(k)'P(k)W(k) + w)ot + 7 


@,(k)P(k I wtk)m(k) +. 7 


~~ (OPCk) Wik) at awim(k) (WC Pio? oo walle awa 


@P(k)@' - SP(k)W(k)m(k)W(k)'P(k)¢' 
- @P(k)W(k) CW(k)'P(k)W(k) + w)-twn' 
- wlwW(k)'P(k)W(k) + w)-'w(k)'P(k)S' 


- nwlw(k)'PC(k)o(k) + w)ltwn' + nwa’ 


@P(k)b' - SP(k)W(k)m(k)W(k)'P(k)®' 
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[-w(k)'P(k)W(k) - w+ W(k)'P(k)W(k) Jn! 
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[W(k)'P(k)wW(k) + w)-tw(k)'P(k)¢' 
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APPENDIX 4B 


Proof of Lemma 1 


The feedback system (4.39) can be described as a 


time-varying autonomous system: 
Ae dae = re CK (ke) (B1) 
where 


came ae 2) oR) 
SR (B2) 
Souk) ease bk) 


ean k) = oF. + ho(k)' Po. (k)e(k) 65 me () (B3) 
Ei2(k) = -h@(k)'m;(k) (B4) 
Eo yt(k) =*-PES (koi ki oetms(k) (B5) 
Sa7tk) =~ 1-1 Pos.(k) otk) Ck) mak) 

=) [T+ Road keuk ihe (B6) 
mice lom(k),pekyeue (B7) 


In equation (B6) we have used the matrix inversion lemma: 
P22(k+1)"' = Po kiero Ck) Ot kaye (B8) 


Now, consider the following quadratic function as a 


candidate for the Lyapunov function: 


Vick )r a=: x CK) TCX CKD (B9) 


where 
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ge Kies 

Tees CK) = T,,(k) = 0 

T22(k) = P22(k)7' 
We obtain 


=a CRAM J) me CK) ay Ce EO CK Ck) 


where Uk were CKoe Gk) = ke) eek) 


(B10) 
(B11) 


(B12) 


(B13) 


(B14) 


A tedious but straightforward matrix manipulation of 2(k) 


yields 
Cah) = SEPELECL I et 
+ (kk)? Poo (k)o(k POs hE. 14+*Fe Lhe dams 
+ [o(k)'Poo(k)o(k) ]26sh Zhesms (ky? 
+ o(k)'P229(k)6262'm3(k) 
es Ck) = -F, Dhd(k)'m;3(k) ad 626(k)'m3(k) 
atl o(k)'P22¢(k)O2h'Zho(k)'m;(k) 
Ope (Uk) =. <o(k) hXSétk ems 0k) = o(k)o(k)'m3(k) 
where 


Qy4(k) 212k) | 
Q(k) = 
Q12(k)' Qe2(k) 


Baia) 


(B16) 


(B17) 


(B18) 
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positive definiteymatrix Zisuch that) (Hitz, iet “al. »v1969): 


Brees bes ee oy ee ES (B19) 
Peon = 50> + rt (B20) 
Deena ese r | (B21) 


Using these relationships on equation (B18), we have: 


M11(k) = -§(kK)S(k)" - O(K)'P22(k)O(k)O262'ms3(k)? (B22) 
Teetk) es > réoCk)? ms,Ck) 
UOC KE ota CR OMEKI Ooo Uke Casares tye ( ee) (B23) 
Geek) so— 9 O UK) OCK) bre eo Cky) oP okey Ok) ) ora ke) (B24) 
where 
fk = £—) o(k) Pes tkieo(k)O57ms, (koe (B25) 
Since 


Q,,(k) a Oe Ck) Qa R Ee Sige Ke 


'P22(k)o(k) 
o(k)'P22(k)¢ {[& + 7r6,m,][$ + 7r02m3]' (B26) 


r? + $(k)'P22(k)¢(k)) 


is semi-negative definite 


and Q2505 22427 = ie (B27) 
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by Lemma 2 in APPENDIX C, the matrix Q2(k) is semi-negative 


definite where "*" denotes psSeudoinverse. It immediately 


follows that the matrix =(k) is stable, so is the system 


Va039):. 
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APPENDIX 4C 


Lemma 2: A given symmetric matrix Q 


1, 242 
Qe = en 
212' N22 ; 


1S semi-negative definite 


ft Q22 < 0 (C2) 
M14 — 24 2%227°212' SF 0 (C3) 
and CeoSiesrQret = OY 2 (c4) 


where "*" denotes pseudoinverse. 


Proof: 


Consider the quadratic form: 


V = Saar) Se un y an me 


= 6 NOT ts 2x'012Y 1 y'Q22y ees) 


where x and y are any non-zero vectors of appropriate 


dimensions. V can be expressed as: 


Ve= pare are < ot 2x D2 (Qee tls 2) y ot y'222Y 
+ x'Q, 2222722 2%22°212'X 


an X12, 2222722222272, 2' 


mex S044) ae 7, 2%22°2;2')x 
PCy, ot Os Gia Bes Ly Qe2* Qe" x) (C6) 


V is non-positive. Hence, 2 is semi-negative definite. 
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APPENDIX 4D 


In this appendix, we introduce an important type of 
stochastic process called supermartingale and discuss those 
of its properties pertinent to the proof presented in 


Section 4.6. 


Definition (Kushner, 1971): Let ix(k)} be a Sequence of 
random variables measurable with paeoeee to. (wert) waerandom 
variable z(k) and {7(k)} a sequence of an ever-increasing 
o-algebras generated by {z(1), --:, z(k)}. The sequence of 


random variable x(k) is then called a supermartingale if 
BUE(k)Ck> ei x (k= 1) 
wacnmprobability 1) (w. pt-i)> for each k. 


One of most successful applications of the 
Supermartingale theory is in the proof of the convergence of 
recursive algorithms for stochastic systems. The 
convergence proof using the supermartingale theory usually 
results in stronger convergence properties compared to that 
in the mean squares sense. The application of the 
eupertatingale theory for the convergence proof largely 


hinges on the following lemmas. 


Lemma 3 (Kushner 1971): If {x(k)} is a nonnegative’ 


Supermartingale sequence and E{x(k)} < @ for all k, then 


Byecck 7 (k-2)3) <x Uk=2) 
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The important implication of the above lemma is that E{x(k)} 
eonverges tova finite random variable. The proof of this 


implication could be found in (Neuve, 1975). 


Lemma 4 (Neuve, 1975): Let {a(k)} and {x(k)} be sequences of 


nonnegative random variables measurable w.r.t $¥(k) and 
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finite random variable x < ©. 


Lemma 5 (Solo, 1979): Let {B8(k)} and {s(k)} be sequences of 


nonnegative random variables measurable w.r.t 7 (k) and let 
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CHAPTER V 
SIMULATIONS 


Daa ENeEToOaduction 

In the previous chapter, three recursive 
identification algorithms for the linear stochastic system 
were derived. The algorithms were obtained by modifying the 
Standard EKF algorithm for identification in order to reduce 
the computational burden and/or to avoid biasedness of the 
estimates. Of the three algorithms, it was shown that the 
third algorithm (EKF-M3) is the most general and hence most 
useful. In this chapter, the effectiveness of this 
algorithm as compared with two other algorithms previously 
reported in the literature will also be discussed with the 


help of simulation studies. 


Since the ultimate objective of developing these 
algorithms is to apply them to real-life systems, it would 
seem appropriate that the comparison of various algorithms 
should be carried out with respect to real-life systems. 
However, this has certain disadvantages. In applying an 
algorithm tomavparticular (system, especially if very Tittle 
oe known about the system, special approaches or tricks have 
to be used. These tricks may vary from algorithm to 
algorithm and system to system. Consequently a proper 
comparison of the different algorithms may become difficult. 


On the other hand, computer simulations provide a great deal 
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8 


of insight into the relative strengths and weaknesses of the 


different algorithms tested for the following reasons. 


i) complete control can be exercised on the conditions 
associated with the system in order to achieve 
clear comparison 

and ii) it is possible to eliminate or reduce dispensable 
elements which may otherwise obscure the comparison 
of results. | 

The comparison of algorithms will therefore be carried out 


by means of computer simulations. 


Many different recursive algorithms have been proposed 
(Saridis, 1974: Isermann, et al., 19742 Soderstrom, et al., 
1978: Dugard, et al., 1980). Listed below are the most 


well-known algorithms which have been published: 


RLS - the recursive least-squares algorithm 

RIV - the recursive instrumental variable algorithm 
RGLS- the generalized least-squares algorithm 

ELS - the extended least-squares algorithm 

RML - the recursive approximate maximum likelihood 


algorithm. 


Considering that both RLS and RIV algorithms are formulated 
to estimate only the process parameters, not the noise 
parameters which are of main interest here, these two 
algorithms will not be used in the simulation studies. 


Although the RGLS method is intended to estimate the noise 
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parameters as well as the process parameters, the model 
Structure representing noise differs from that used in the 
EKF-Maiiand «ther last two valgorithms listed. ‘Consequently, 
the RGLS method will be of little value in the comparison 
Studies in terms of assessement for possible real 
applications; thus this algorithm is also excluded. The ELS 
and RML algorithms will be used in simulations along with 
the EKF-M3 algorithm developed in this thesis. It is noted 
that the features of the ELS and RML algorithms were 


dicussed in the previous chaper. 


Two examples will be used for the simulation studies. 
In Section 5.2, results based on a well known example which 
first appeared inj (Ljung, et al., 1975) are discussed. 
These studies will compare the convergence properties of the 
algorithms. In Section 5.3, general comparison of the 
algorithms is reported based on a fourth order example used 


extenSively by Lee (1964) and Saridis (1974). 


5.2 First Simulation Study 


A second order system is chosen: 


ay 1 da; 
x(k+1) = pa) bce v(k) 
a2 0 d2 ; 


y (Ry = 1 [1 Od Qik Bai) 
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and -v(k) ~ N(0,1) 


This example is the same as that used by Ljung, et al. 
(1975) as a counterexample to demonstrate a biased 
convergence of the ELS algorithm. Since the EKF-M3 
algorithm as well as the RML algorithms yield unbiased 
estimates, the ELS Aue tie was modified also to yield an 
unbiased estimate so that a true comparison of rates of 
convergence of the three algorithms could be made. The 
modified ELS algorithm is designated as GELS. The details 


of the modifications are given in Appendix 5A. 


Several values of initial state were used in the 
Simulation studies. While the GELS and the EKF-M3 algorithm 
converged every time with a proper choice of vector f, the 
convergence of the RML algorithm was possible only by the 
useso£— a projection facility such that estimate c(k) stays 
in a region where 


1 


1 2.e0(k jase -fchbk)q-4 


is always stable. 


Finally for the EKF-M3 and GELS algorithms, a vector f 


was selected as: 
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Figures 5-1 and 2 show the results as well as the 


measurement y(k). 


229 Second Simulation study 
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For the second comparison study, the following fourth 


System is chosen: 


k+ 1) = Ax@h)asaro Cio) 


k= A x CK), pies) 


— ai 


h C1 OTRO Ons 


the Gaussian noises 
wk) ~ NL. OF oO) 


ytk) ~sNCO0L Oar) 


For purpose of simulation, the parameters to be identified 


abe set ass 


This 


— 
e 
©o 
(5) (a(S) (2) 
e e e 
(2) 
oO 
— 
e 
Oo 
Ww 


00° = 182 


i) 
Oo 
N 
Ol 


system is equivalent 


opel 


OGrlst 


. 60 
.00 


89 


See 


terms of input-output relation 


85 


sd ae for oe oni 


iz2wo0) posse lloetam . geute y ebaac 


ea . ' . ’ ee i e 


my ba io oy bs Bri 
78g | a Dred sai 

batt amas adh ag Pacem: ae 

ieee ie a 


2 atta Mie 7% eae 
ae ata aaaaa 1s sai 


aotseley JuAseo*tugEl To WM Be soolevinpe ai wsdaye ke eid 


thet 
ee 


°ae bite ai 


y 


86 


0002 


(¥)4 Juaswainseaw ayy JO una 4nd4yno sayuL 1°s 


djls 
00st o00t 00s 


00° OT- 


00°0 
(MJA LNSWSYNSBUSAW SHL 


00° OT 


- rs ; 
Lnmy O€ 


onciBinr 


; Lv 


5 


- 


‘ ; i 4h, 
i =>) i oe 


em C7 | od Ee 


fas fF 


a 


o-00 CO 9°00 


WE WEBRNYEWEML ALK) | 


87 


FD) JO Sezer jso s4Lr ec SG 


cae 
0002 Oost o00t 00s 0 


00° 2- 


00° T- os*T- 
SLBWILSS 


os°0- 


00°0 


iv 
er 
ros . 
— 
a 
oi 


>a 


2 


88 


7) JO saqeutqsa auL qz°s 


aac S 
0002 00st ooo 00s 0 


O92 1= 


00° T- 


00°0 0s°0- 
JLBWILSS 


os°0 


89 


'@ JO Sazeutjss ay, oz‘c 


djls 
0002 00st o00t 00s 0 


OS: T= 


00° T- 


00°0 os*0- 
SLBWILSS 


os°0 


i 


— 
Sat rugres 


upe 


2 '.3- 


90 


te Jo sazewtqss syuL P7’s 


cena 
0002 Oost ooot 009 0 


00°0 abticl thes 
3LUWILSA 


os°0 


2S © 
ad i 


were 


Ser 


Sy Efe 


~~ 
i 


ca 


ll 
= 
= 


oleb 


PA Sn) Pa eaw e'oh G'S) me Liars hia Ge) 


y(k) Nee CK ie teu k') 


mere C1 =1 00,570 =0.027 0306 =ORo7ionhe 
thus pres—r[-044304— 0153 Qed Tdoe— 0-0 4ete}e 


and v(k) ~ N(0.0, 3.988) 


Two featuresmetethis system are worth noting. This 
System is very sensitive to noise (Lee, 1964) and 
Gonsequently, .tawasedutiacultetowidentify the system 


parameters without bias (Saridis, 1974). Also 
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The same three algorithms used in the first simulation 
Study are also used here. For the EKF-M3 and GELS 
algorithms, the following values were assumed as the 
a@ priori informations on the vector parameter c. The same 


values are also used as the vector f: 
CO) Sect eee 210470--20.52.0...3.0 0.30) 


The RML algorithm was also initiated with the same initial 
values for fair comparisons. The RML algorithm is also 
complemented with the projection facility as the previous 


simulation study. The results of the simulation studies are 


shown in Figure 5.3 and 5.4. 
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5.4 Remarks 

The results of the comparison are summarized in 
Rigures 5-9 and 5e6 for both studies.. These figures 
illustrate the convergence properties of the algorithms. 
The squared estimation errors were normalized by dividing 
then by the initial squared errors to give the same weight 


to each estimate. 


Both studies show that the performance of the EKF-M3 
was vastly superior to those of the other two algorithms for 
almost all the the parameters. The performance of the RML 
algorithm on the other hand varied considerably depending on 
the parameters. For example, in the second study its 
convergence. rates for estimates c,(k) and c3(k) were 
unacceptably slow while those for c2{k) and cy(k) were 
Surprisingly fast. Though the overall convergence rate of 
the RML algorithm was not much slower than other two 
algorithms as shown in the figures, this algorithm 
experienced the worst transient reponse among the three. As 
a result, even after 5000 iterations the convergence of this 
algorithm was not conclusive. As expected, the GELS 
algorithm showed slower convergence than the EKF-M3 
algorithm. However the convergence patterns of the both 
algorithms were quite similar as shown in Figure 5-4 because 
the gradients of the error correction terms of the both 


algorithms are the same. 
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APPENDIX 5A 


Derivation of the Generalized Extended Least Squares 
Algorithm (GELS) 

AS was pointed out in Section 4.5, it should be noted 
that the difference between the ELS algorithm and EKF-M2 
ecoocithm 1s /thateun ELS a priori’ errombe tk), 1s: used tor the 
generation of 6,(k) while the EKF-M2 algorithm uses 
aepesteritorirtérror eiCkiSecl feshoularalso maei 62d that the 
EKF-M3 algorithm was obtained by generalizing EKF-M2 using a 
canonical form which included a general vector ftito avolide 
the possibility of biased convergence. These two ideas will 
now be used to obtain a generalized Extended Least Squares 


algorithm. The GELS algorithm is then of the form: 


zo(k+1) = F'z2(k) + hy(k) (Aa 
6,(k+1) = P'6,(k) + hly(k) - $(k)] i (A1.2) 
D(k+1) = Blk) + Po2(k)o(k)m3(k)[y(k) - ¥(k)] GAT AS: 
y(k) = ¢(k)'p(k) GA Aa) 
o(k) = [6,(k)",22(k)"]! (A1.5) 
Poo(kt+1) = Poz(k) - P22(k)o(k)ms(k)o(k)'P22(k) (A1.6) 
Bem isonet poi ke UD oer ake ont) a (AE7) 


P,,(0) = P22(0)' eG 


Covergence of GELS 
The errors p(k) = p - p(k) and 6,(k) = 6,(k) - 64(k) 


satisfy the following error equations: 
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6,(k+1) = FYO,(k) - hy(k) CAD.) 
e(k) = 62'6,(k) + y(k) + v(k) (22) 
Doki = DCk)e = Po (kok) ms, (Kye OR) CA 23) 
y(k) = ¢(k)'p(k) (A2.4) 
Poke n= P22(k) -r Po o(k)eUk)ms (Ck) ok)! Pook) (A2y 5) 


Whevenro. =..F —' Osh! 


As was the case for the EKF-M3 algorithm, the GELS algorithm 
is also assured of unbiased convergence if the errors p(k) 
and @,(k) converge to the origin in the error vector space. 
Furthermore, if the errors p(k) and 6, (k) converge to the 
origin without the explicit noise v(k) in the equation 
(A2.2), one can easily conclude by similarity of the 
structure of the error equations that the GELS algorithm 


will yield unbiased estimates of the parameters. 


Consequently, the overall unbiased convergence of the 
GELS algorithm can be established by examining the stability 
of the error feedback system without the noise term v(k) in 


equation (A2.2)(«)Thissstabilaty  yproblem@cantbeystated as 


follows: 


Lemma: Given the following time-varying system: 


G8, (k+1) = FYO,(k) - hy(k) (Asus 
e(k) = @2°6,(k) + y(k) (A3.2) 
D(k+1) = p(k) - P22(k)o(k)m3(k)e(k) (A3.3) 


y(k) = (k)'plk) | (A3.4) 
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Po2(k+1) = P22(k) - Po2(k)o(k)m(k)o(k)'P22(k) 


P22(0) = P22(0)' au 


6,(k) and B(k) converge to zero as k > @ 
PreG(kK) Peek) @Ck) |S 5 ‘for all k where 6 15-4. positive 


Scalar and 
Huan wt Aoi Go oe 1 SiR e pb) ib 
is strictly positive real. 


Proof: 
Consider the following quadratic function as a 


candidate for the Lyapunov function: 


Vik) Vike Ve woke) 


~— 


Vack) = 6,(k)'r6,(k) 


V5.0k)) = DURY Poet kee pk) 


It follows from equations (A3.1-2): 


Walk ess Viel kT) eek) 
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VVi5.k) = ViosGket a) a V>.Gk) 


y(k)? - 2y(k)e(k) + @(k)'P22(k)o(k)m3(k)e 


In equation (A8), the following matrix inversion lemma 


been used: 


(A3.5) 
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Sce Hige NG =a 2's 6oial = Faye thmismsap yon, there sextet 
ea realsvector ¢, ai real scalar r and a real positive 
definite) matrix 2 such, that: (Hitz, .et al..),1969) 

FLIF,' =D = -&¢! (A10) 

|carap oN gh = 05 sy tT CASIhIe) 

Mies hl ee (A12) 
With these relationships, equation (A7) becomes: 

Waki) 6 er Pm Le eek Bite Ck) or itt ete ea lncuyy wiry BI )iee ANS) 
thus, 

VV (k) = VV—-Uk) as VV2(k) 

= -{[6,(k)'& + y(k)7TJ]? - $(k)'P22(k)o(k)m3(k)e(k)?} (A14) 
Since ¢(k) is bound and P22(k) is nonincreasing, there 
exists a positive scalar 6 such that 

$(k)'Po2(k)o(k) < 6 (A15) 
in which case VV(k) <= 0, 


and the feedback system 


(A3) 


is stable. 


Consequently, the GELS algorithm is assured of producing 


unbiased estimates of the parameters. 


The slower 


convergence of GELS compared to EKF-M3 can be attributed to 


this additional restriction (A15). 
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CHAPTER VI 
CONCLUSIONS 


6.1 Summary 

This thesis is entirely devoted to the development of 
Various recursive identification methods for deterministic 
as well as stochastic discrete linear time-invariant 
Systems. Though at the beginning of each chapter, there was 
a brief discussion on the specific subject which was dealt 
in that chapter, it is the purpose of this section to 
Summarize the salient features of each chapter and the 


results obtained in this thesis. 


In Chapter II, a canonical form of the innovations 
representation for the multiinput singleoutput system is 
introduced. The preference for using the innovations 
representation over other representations, in particular for 
the identification problem is largely related to the 
uniqueness of parameter estimates for a given measurement. 
This uniqueness, which is of great importance for the 
Pacee er eas problem, originates from the fact that the 
innovations representation is the sole minimum phase system 
for a given rational spectral density (Anderson, et al., 
1979). This, in turn, implies that for a given measurement 
y(k), a stationary process, there exists only one minimum 
phase system whose output could be regarded as the 


measurement and this minimum phase system is identical to 
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the innovations representation. Thus, by adoption of the 
innovations representation, a risk of ambiguities of 
identified parameters could be almost eliminated. 
Furthermore, the structure of the canonical form introduced 
is Suited especially to the least squares method for the 
identification problem. This form is also flexible enough 
to facilitate the development of new identification 


algorithms. 


Chapter III presents a fast identification algorithm 
for the deterministic linear system. This was done by 
geometrically interpreting the Kudva-Narendra (K-N) 
identification algorithm and improving the convergence speed 
by the introduction of orthonormal vectors. The special 
advantage of this scheme is that the overall computational 
requirements are considerably reduced in comparison to the 
K-N algorithm. While the underlying idea of this 
development is simple and straightforward, the geometric 
interpretation makes it possible to visualize how the K-N 
algorithm works and consequently to improve the convergence 
rate. Furthermore, this approach provided a basis for the 
direction to be taken for the development of identification 
algorithms when noise is present, i.e. the stochastic system 


identification which is the subject of Chapter IV. 


A stumbling block in the use of the standard least 
Squares method for the stochastic system identification 


problem is, in general, the inclination of the least Squares 
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method to bias the estimates (Astrom, et al., 1971) if noise 


is correlated, as is generally the case. Several techniques 


have been proposed for the recursive identification problem 


to get around this drawback. The basic philosophy behind 


the proposed methods in common is to decompose the 


correlated noise into a sequence of weighted white noise or 


innovations, which is often called as Wold decomposition and 


then to use the standard recursive least squares method to 


identify the augumented parameters which includes the 


parameters for weighing the innovations. One of the obvious 


problems arising in this approach is that the innovations 


are not amenable to meaSurement. Consequently, different 


methods have been suggested to estimate the innovations. 


Since the identification problem for the linear 


stochastic system can be considered as a nonlinear 


estimation problem, the Extended Kalman Filter algorithm has 


been quite often used as a starting step for developing new 


algorithms. In chapter IV, this approach is once again 


taken to derive several different identification algorithms 


in order mainly to alleviate the unnecessary computational 


burden. The modifications are made possible because of the 


following observations: 


i) 


2) 


the EKF algorithm for the identification problem is 
assured of its unbiased convergence in the neighborhood 
of the true values of the parameters. 


the use of the innovations representation combined with 
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the above observation justifies setting some elements of 
the approximate covariance matrix of the EKF algorithm to 
be zero, thus leading to Simpler algorithms in term of 
computational requirements. 
The process of simplification has also revealed that some 
existing algorithms which were developed by other methods 
are merely modifications of the EKF algorithm. Thus, this 
helps to understand better the behaviour of these 
algorithms. This understanding together with an effective 
use of the canonical form introduced in Chapter II makes it 
possible to develop a fast but bias-free identification 
algorithm if some condition is met. The requirement of the 
condition could be in most cases satisfied with a priori 


knowledge of the system. 


The differences among the algorithms discussed in this 
thesis, namely, RML, GELS and EKF-M3 are worth noting. 
While in the RML and GELS algorithms, a priori errors are 
used to replace the true innovations, a posteriori errors 
are used for the EKF-M3 algorithm. It is this a posteriori 
error that is thought to be responsible for faster 
convergence of the EKF-M3 algorithm to the true values 
compared to the RML and GELS algorithms. Secondly, the RML 
algorithm incorporates a prefilter with time-varying 
parameters. On the other hand, the GELS and EKF-M3 
algorithms employ a prefilter with fixed parameters. 


Theoretically, identification algorithms with time-varying 
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prefilter exhibit superior performance compared to 
algorithms with time-invariant filters, provided the 
parameters of the time-varying prefilter approach to the 
final desired values uniformly. In the case of the RML 
algorithm this uniformity property is not present in the 
prefilter. Consequently its convergence property is 
inferior to that of either the GELS or EKF-M3 algorithm as 
shown by the simulation studies in this thesis. Table 6-1 
summarizes the properties of the various recursive 


algorithms discussed in this thesis. 


It should also mentioned that the analysis of Ljung on 
the convergence of the RML algorithm (Ljung, 1979) though 
elegant, proves just the mere convergence of the algorithm 
and does not address the question whether the convergence is 
either uniform or asymptotic. One of the advantages of 
using a time-ivariant prefilter as in the EKF-M3 algorithm, 
is that it makes the analysis of the convergence property of 


the algorithm straightforward. 


Though the most of discussions have been for 
multiinput singleoutput systems, the extension for 
multiinput multioutput systems can be easily made (Anderson, 


1974; 1977). 
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6.2 Suggestions for Further Research 

Fast but bias-free identification algorithms for the 
linear system using a systematic approach have been 
developed. But the ever-increasing use of on-line 
identification methods for systems in the various fields 


substantially expands the range of application of the 


identification techniques to systems in which more stringent 


conditions are often encountered; for example, it is often 
necessary to identify a system with very limited a priori 

knowledge of the system. Thus there is a need to improve 

the existing algorithms so as to handle this situation 


Satisfactorily. 


The employment of a time-varying prefilter instead of 
a time-invariant one used for the algorithm developed in 
this thesis seems extremely promising as an immediate 
improvement, thus eliminating the requirement of a priori 
information to guarantee unbiased convergence. One very 
important problem requiring an immediate attention is the 
unbiased convergence property of the modified form of the 
algorithm. This problem also demands a comprehensive study 
of a more general form of representation for the stochastic 


system to accomodate the time-varing prefilter. 


= 


Another approach which has a great potential is to 
treat the linear stochastic system as an unknown infinite 


impulse response (IIR) filter and use the standard least 
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Squares method. This approach will involve the development 
of an efficient recursive least squares method for an 
infinite dimensional system and a fast algorithm to 
transform a polynomial of ever-inceasing order into a 


rational polynomial with a fixed order. 


While this thesis has been devoted to the theoretical 
development of efficient recursive identification algorithms 
for stochastic systems, it 1S appropriate to mention that 
these algorithms have a wide application in the area of 
process control (Isermann, ed; 1980), communications 
(Friedlander, 1982) and processing of seismic data 
(Robinson, «957: Mendel. 1977>Mahalanabiyss etial.. 1981). 


obtained from geophysical exploration to mention a few. 


More extensive computer simulations also needed for a 
more complete evaluation of the algorithms treated in this 
thesis. A comprehensive assessement of the performance 
characteristics including sensitivity to error in estimation 
system order and the effect of unusually large impulse 
noise, etc., which are of extreme importance in real 
applications should be made by means of simulation studies 


before implementation. 
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COMPUTER PROGRAMS 


Following are the listings of the computer programs 
used in Simulation studies of Chapter V. Since, as 
mentioned previously, the EKF-M3 and GELS algorithms differ 
only in the use of a postériori or a priori errors, these 
two algorithms are combined together into one and selected 
by a flag ICH of the routine GELS. The calculation of the 
gain vector of the correction term of the parameter 
estimation is performed with the fast algorithm proposed by 
Ljung, et al. (1978). The program GELS includs this fast 
algorithm. The rest of GELS program and entire RML program 


follow exactly the same notations as used in Chapter IV. 
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SUBROUTINE GELS(YU,FV,P,ITER,1YU,NORD,NORDD,GAIN,ICH) 


THESE ARE THE GENERALIZED EXTENDED LEAST SQUARES AND 
EKF-M3 IDENTIFICATION ALGORITHMS FOR THE SYSTEM 
WITHOUT CONTROL INPUT, i.e. TIME-SERIES, OR WITH 
CONTROL INPUT. 

THE CONTROL FLAG 'ICH' SWITCHES BETWEEN GELS AND 

EKF 4M3% 

THE RESULTS WILL BE STORED IN A TEMPORARY FILE "-GELS" 
Pe .GELSSLS (RUN TORG) SEKEoM3 JOR ITBKE-M3i°1S’ RUN. 


YU; OUTPUT OF THE SYSTEM WHOSE PAPRAMETERS ARE TO BE 
IDENTIFIED AND THEN THE CONTROL INPUTS ARE 
IMMEDIATELY FOLLOWED IF THEY EXIST. 


FV: SYSTEM VECTOR OF THE SYSTEM MATRIX F OF THE 
IDENTIFIER 
T 
I 
Fi = BV tiers 
0 
P: THIS VECTOR CONTAINS THE INITIAL ESTIMATES FOR 


¢, a and b VECTORS 
ITER: THE NUMBER OF ITERARTION DESIRED. 
EU's THE NUMBER OF ENTRIES IN "YU" IF CONTROL INPUTS 
: EXISTS THEN IYU=ITER*2, OTHERWISE IYU=ITER 
NORD: ORDER OF THE SYSTEM 
NORDD: 2*NORD IF NO INPUTS 
3*NORD IF INPUTS EXIST 
GAIN: SCALAR VALUE FOR THE INITIAL GAIN MATRIX P22 
P22(0) WILL BE GAIN*IDENTITY MATRIX 
ECHs SELECT GELS IF ICH=0 
EKF-M3 IF ICH=1 


POGI CALS) DATA 1 (6) DATA2 (39 

READ EK (1S) (115) SACS, 8): EVINORDIG Digi o53) PX IG3)s 
PRMT( 15), YUCLYU) ,EPST.(3)))WKAG3, 4) 1STG(353)5 
AMAX(15),AMIN(15),P(NORDD) ,EPSI0(3) 


DATAGDATANKL= "246Gb e" LS! oe / 
DATA DATA eae ks py ee oa, by 
ATTACH THE PROPER TEMPORARY FILE FOR THE OUTPUTS 


IF(ICH .EQ. 0)CALL FTNCMD('ASSIGN 3=?;',0,DATA1) 
IF(ICH .EQ. 1)CALL FTNCMD('ASSIGN 3=?;',0,DATA2) 


DETERMINE IF THERE ARE INPUTS TO THE SYSTEM 
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NDATA=NORDD/NORD 


SET SIG(0)=GAIN*IDENTITY MATRIX. THIS IS EQUIVALENT TO 
SETTING P22(0)=GAIN*IDENTITY MATRIX. THUS, THE INVERSE 
OF SI1G(0)), -I1SIG(0)=IDENTITY MATRIX/GAIN 


DO 1 I=1,NDATA 

DO -2.J=1,NDATA 

ISIG(1,J)=0.0 

BRAGRMCEO STI BISIG(1, J) =1.0/GAIN 
CONTINUE 

CONTINUE 


SED T/A €0)| ANDED(0) EQUAL TO ZERO 


DO 5 I=1,NORDD 
DO 6 J=1,NDATA 
A Glen y= O00 
BGI IH 6.0 
CONTINUE 


SET (0). 7ANDexiCO)) EQUAL: TO ZERO 


SET INITIAL ESTIMATES P, P=[FV-C,A-FV,B] 


DO 8 I=1,NORD 
P(I)=FV(I)-P(I) 
P(I+NORD)=P(I+NORD)-FV(I) 


SET INITIAL VALUE OF AMAX AND AMIN EQUAL TO ZERO 


DO 9 I=1,NORDD 
AMAX(1)=-100.0 
AMIN(I)= 100.0 


DO 1000 KKK=1,ITER 


YH=0.0 

L=1 

DO 10 I=1,NDATA 

DO 11 J=1,NORD 
YH=YH+P(L)*X((J-1)*NDATA+I ) 
L=L+1 

CONTINUE 
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ERR=YU (KKK) -YH 
SEEebgn. (4.13) in) Ljung's Paper in TJC 


AH=1.0 
bee. CHI ¥EOLLOJGOTO 12 
AM3=0.0 
DO 15 I=1,NORDD 

15 AM3=AM3+K(I)#*X(I) 
AH=AH-AM3 


t2 DON T78 IS 1., NORD 
PRMT(I)=FV(I)-P(I) 
PRMT(I+NORD)=FV(I)+P(I+NORD) 
iE ACNDAITAY, © BO... <2) GOTO. 17 
PRMT (I +NORD*2) =P(I+NORD#*2) 
17. CONTINUE 


DO 18 I=1,NORDD 
AMAX(I )=AMAX1(AMAX(1),PRMT(I) ) 
18 AMIN(I)=AMIN1(AMIN(I),PRMT(I) ) 


WRITE(3,200) (PRMT(J) ,J=1,NORDD) 


200 FORMAT(1X,15F12.4) 


L=1 
DO 20 I=1,NDATA 
DO 21 J=1,NORD 
P(L)=P(L)+ERR#*K((J-1)*NDATA#+I ) 
24 Gebt4 
20 CONTINUE 


CALCULATION OF NEW COMER FOR SHIFT REGISTER "X" 


DO 25 I=1,NDATA 
S=0.0 
DO 26 J=1,NORD 
26 S=S+FV(J)*X((J-1)*NDATAHI ) 
-IF (I .EQ. 1)HIN=ERR*AH 
tr {1 JEO.N2)HIN=YU(KKR) 
IF (I .EQ. 3)HIN=YU(ITER+KKK) 
25 STAT) =S+HIN 


KHKKKKKKKKKKKKKKKKKKKKKKKKKK KKK KK KKK KKK KKK KKK KKK KKKKEKKEKKK 


* THE ITERATIVE METHOD OF CALCULATION OF "K" * 
* “BY E5 ung) TSMUSED. SEE Int. J. of Contr. Vol. 27 No 1+ 
* -ppsiaid,v 197s" * 


KKK KK KK KK KKK KR KKK KK KK KKK KK KKK KK KK KKK KKK KK KKK KK KK KKK KKK KK 
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CALCULATION’ OF’ "EPSILONO”” OF Eqn.” 4.11 


DO 30 I=1,NDATA 

S=0.0 

DO 31 J=1,NORDD 

S=S+A(J,1)*X(J) 

EPSIO(1I)=XI(1)+S 


URDATEe As (OF eqn.) 4:.)2 


DO 33 I=1,NORDD 

DOo 34i J=1,NDATA 
A(i,d)=A(1,0)=K(1) #EPSI0(J) 
CONTINUE 


CALCULATION OF "EPSILON" OF Eqn. 4.3 


DO! 40=1=1*)NDATA 
S=0.0 

DO 41 J=1,NORDD 
S=StA(J,1)*X(J) 
EPSI(1I)=XI(1I)+S 


UPDATE "ISIG" OF Eqn. 4.4. HERE THE MATRIZ INVERSION 
LEMMA HAS BEEN USED. 


S=1.0 

DO 90 I=1,NDATA 
WKA(I,1)=0.0 
WKA(I,2)=0.0 


DO 91 J=1,NDATA 
WKA(I,1)=WKA(I,1)+ISIG(I,J)*EPSI (J) 
WKA(I,2)=WKA(I,2)+EPSIO(J)*ISIG(J,1) 
S=S+WKA(I,2)*EPSI (1) 


DO 95 I=1,NDATA 
DO 96 J=1,NDATA 

PSIG(IL 9) =—SIGUL, J) -WKA( 1, 1) 4WKA(I,2)/7S 
CONTINUE 


CALCULATIONY OPS K>bar/ OF -Raqnse 475 


DO 45 I=1,NDATA 
S=0.0 

DO 46 J=1,NDATA 
S=S+ISIG(I,J)*EPSI (J) 
WKA(I,1)=S 


DO 48 I=1,NORDD 
S=0.0 

DO 49 J=1,NDATA 
S=StA(1I,J)*WKA(J,1) 
K(I)=K(1I)+S 
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DK PAS te M(t sine Ronee 6 * 


KKK KKK KKK KK KKK KKK KKK AK KKK KKK KKK KKK KKK KKKEKEKKEKE 


50 


a 


52 


56 
Shs) 


58 


61 
60 


66 
65 


UPDATEmex” AND COMPUTE) "M(t)" and) "MU(t)" in Ban. 4.6 


NN=NORDD-NDATA 
DO 50: 1=1;/NDATA 

EPSI (I )=K(NN+I ) 
WKA(I,2)=X(NN+I ) 

DO 51 I=1,NN 

K (NORDD-I+1)=K(NN-I+1) 
X (NORDD-I+1)=X(NN-I+1) 
DO 52 I=1,NDATA 

X01) =X1 1) 
K(I)=WKA(I,1) 

XI (1) =WKA(I , 2) 


CALCULATION OF "BTAOQ”" OF Eqn. 4.7 


DO 55 I=1,NDATA 

S=0.0 

DO 56 J=1,NORDD 

S=S4D(T plex (S) 

EPSI0(1)=XI-(1)+S 


URBATES"D® OF Eqn. 4.8 
FIRSTLY, THE SECOND TERM OF THE RIGHT SIDE IN Egn. 
AGAIN THE MATRIX INVERSION LEMMA HAS BEEN USED. 


S=1.0 
DO 58 I=1,NDATA 
S=S-EPSI (I) *EPSIO(I) 


DO 60 I=1,NDATA 
DO 61 J=1,NDATA 
WKA(I,J)=EPSI(1I)*EPSIO(J)/S 

IF (I .EQ. J)WKA(I;d)=WKA(1,J)+1.20 
CONTINUE 

CONTINUE 


THEN, THE FIRST TERM OF THE RIGHT SIDE IN Eqn. 4.8 


DO 65 I=1,NORDD 
DO 66 J=1,NDATA 
D(I,J)=D(1,J)-K(1I) #EPSI0(J) 
CONTINUE 


COMBINE THE ABOVE TWO TERMS 
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DO 70 L=1,NORDD 
DO 71 I=1,NDATA 


WKA(I ,4)=0.0 


DO 72 J=1,NDATA 


WKA(I ,4)=WKA(1I,4)+D(L,J)*WKA(J,1) 


CONTINUE 


DO” 73 M=1,NDATA 
D(L,M)=WKA(M, 4) 


CONTINUE 


UPDATE "K" 


DO 75 I=1,NORDD 


S=0.0 


DO 176 J=CINDATA 
S=S+D(I,J) *EPSI (0) 


K(I)=K(I)-S 


CONTINUE 


WRITE(3,200) (AMAX(J),J 
WRITE(3,200) (AMIN(J) ,J=1,NORDD) 


WRITE(3,300) 


FORMAT(1X, "ESTIMATE OF C VECTOR IS FOLLOWED BY THAT 


SHOFICAY) 
WRITE(3,301) 


1,NORDD) 


FORMAT(1X,'IN FORMAT F12.4') 


WRITE(3,302) 
FORMAT(' ') 


IF’ (ICH MEO. 
WRITE(3,304) 


FORMAT(1X,'THESE ARE OUTPUT OF EKF-M3' ) 


RETURN 
WRITE(3,305) 


FORMAT(1X,'THESE ARE OUTPUT OF GELS' ) 


RETURN 
END 


0)GOTO 303 


128 


rp 


a 


rt 
an. 7 


Ta ra? 


AOA AOI OGA OAM OA OAM OAAMWA AO) 60'O.0°O 


OO a QXOLCa@a 


YOO 


$ 
$ 
$ 


ic 


SUBROUTINE RML(YU,FV,P,ITER,IYU,NORD,NORDD,GAIN,ICH) 


THIS IS THE RECURSIVE MAXIMUM LIKELIHOOD 
IDENTIFICATION ALGORITHM COMPLIMENTED WITH 

PROJECTOR WHICH PROJECTS THE C ESTIMATES INTO THE 
PREDETERMINED AREA FOR THE SYSTEM WITHOUT CONTROL 
INPUT, i.e. TIME-SERIES OR WITH CONTROL INPUT. 

THE RESULTS WILL BE STORED IN A TEMPORARY FILE "-RML" 


Vai: OUTPUT OF THE SYSTEM WHOSE PAPRAMETERS ARE TO 
BE IDENTIFIED AND THEN THE CONTROL INPUTS ARE 
IMMEDIATELY FOLLOWED IF THEY EXIST. 


EV: THE VECTOR FV IS ALL ZERO FOR RML, THUS FORCED 
ROBE LERO. j 
PS THIS VECTOR CONTAINS THE INITIAL ESTIMATES FOR 


Gy aeandeb VECTORS 

ire Rs THE NUMBER OF ITERARTIONS REQUIRED. i.e. THE 
LENGTH OF ARRAY OF Y 

WEY: s THE, NUMBER: OF ENTRIES IN "YU" IF CONTROL 
INPUTS EXIST THEN IYU=ITER#*2, 

OTHERWISE IYU=ITER 

NORD: ORDER OF THE SYSTEM 

NORDD: 2*NORD IF NO INPUTS 
3*NORD IF THERE ARE INPUTS 

GAIN: SCALAR VALUE FOR THE INITIAL GAIN MATRIX P22 
P22(0) WILL BE EQUAL TO GAIN*IDENTITY MATRIX 

veHs ICH IS ADDED FOR THIS ROUTINE TO BE FORMAT 
COMPATIBLE WITH GELS © 


LOGECALS INDATA 148) 

REARTEO(5e5) }H(5),P22( 15715) [UP G5..15)s, UN 5115) 5 
ZACS 2265) 4231-5) PRMD So hy (io) coal Sack (oom, 
YU(IYU) , TMAX(15),TMIN(15),FV(NORD) ,P(NORDD) , 
CINETUS) BFCC5) 5) 

DATA H/1.0,4#0.0/ 

DATA UP/75*0.0/ 

DARRMOATATY 2-— ly UR YIM! 4c Lye hack ie: 


ATTACH THE PROPER TEMPORARY FILE FOR THE OUTPUTS 
CALL FTNCMD('ASSIGN 3=?;',0,DATA1) 


DETERMINE IF THERE ARE INPUTS TO THE SYSTEM 


INPUT=0 
IF (NORDD .EQ. 3*NORD) INPUT=1 


SET UP FO 


DO 600 I=1,NORD 
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DO 610 J=1,NORD 


BO City Jp= 0.20 

DEG wEOL oer POCT 3) =1.0 
CONTINUE 

CONTINUE 


DO 620 I=1,NORD 
FV(1)=0.0 


SET INITIAL GAIN MATRIX P22(0)=GAIN*IDENTITY MATRIX 


DO 1 I=1,NORDD 
DOwZ) an, NORDD 
P22( 07s) = 0.0 


TPrClg. BOs Je P22( 140) =GAIN 


CONTINUE 
CONTINUE 
DO 5 I=1,NORD 
W341 )=020 
Z2(1)=0.0 
Zeer) =O 0 


SET INITOAL GST MATES? Py, 


DO 7 I=1,NORD 
P(I)=FV(1)-P(1) 


P(I+NORD)=P(I+NORD)-FV(I) 


P=[FV-C,A-FV,B] 


STORE THE INITIAL VALUES FOR ESTIMATE C 


DO 8 I=1,NORD 
CINTT(CT)=P(1) 


DO 650 I=1,NORDD 
TMAXIOT )=—1100.0 
TMIN (1) =30100.0 


DO 77 (TK= IWaTER 


DO 10 I=1,NORD 
BY (1) = 241 (Ty) 
FY(I+NORD)=Z2(I1) 

TF) (INPUT “EQ. 0) GOTO 
FY (I+NORD#*2)=23(I1) 
CONTINUE 
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YH=05..0 
DO 20) I'=1),NORDD 
20 YH=YH+FY(I)*P(I) 


ERR=YU(K)-YH 
CALCULATION OF FC = FO - H*HTHETA2 


DO 40 I=1,NORD 

DO 7420tu.=7) , NORD 
42. FC Bi Jer 0 @0) =H (1) +P (9) 
40 CONTINUE 


CALCULATION OF UP(k+1) in eqn. (4.24) 


DO: 410 T=1,NORD 

DO 420 J=1,NORDD 

UNC) =07%0 

DO. 430) L=1,NORD 
230 UNC H=UN Geet FC (1 ,L) *UPAL,J) 
220 UN GL AR=UN Gl D-H )ieey (J) 
410 CONTINUE 


CALCULATION OF S = P22(k)UP(kt+1)'h in eqn. (4.28d) 


DO 450 I=1,NORDD 

S(1)=0.0 

DO 460 J=1,NORDD 
460 S(I)=S(1)+P22(1,J)*UN(1,J) 
459 CONTINUE 


CALCULATION (Ofdimi(k) with w=1 "in eqn. (4.27) 


AM1=1.0 
DO 58 I=1,NORDD 

58 AM1=AM1+UN(1,1)*S(I) 
AM1=1.0/AM1 


UPDATE 21(k) 


DO 70 I=1,NORD 
TT) =0 30 
DO 72 J=1,NORD 
72 MENSrG EO Ceo Ze (os) 
70) TCM =T(T +H (1) *ERR 


DO 74 I=1,NORD 
TAe ODA) Tle) 


UPDATE 22(k) 


DO 80 I=1,NORD 
THAW Y= 0.10 
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DO 82 J=1,NORD 
RCI em (yar Otte) *Z22(7) 
TH IDET +1) *YU(K) 


DO 84 I=1,NORD 
22.0 yenCls) 


UPDATE 2Z3(k) 


DO 86 I=1,NORD 

Per) =0.40 

DO 87 J=1,NORD 
TOAORaT CIDER OCTED) +23.(99 
Tea) =2 (1) +H 1 ) *YU(K+1 TER) 


DO 88 I=1,NORD 
Zoey) eT (7) 


UPDATE THE ESTIMATE P(k) in eqns (4.28a) and (4.28d) 


DO 90 I=1,NORDD 
P(I)=P(1I)-S(I)*AM1*ERR 


TEST? FORIJ ESTIMATE FOR VECTOR C 


CALL RAIBLE(P,NORD,ITEST) 
ie reST..-EO. 0)GOTO 102 
DO 110 I=1,NORD 

Ph) =CENL? 1, 


DO 550 I=1,NORD 
PRMT(1I)=FV(I)-P(I) 
PRMT'(I+NORD)=FV(I)+P(I+NORD) 
IF PCINPUT TEQ« 0) GOTO 550 
PRMT (I +NORD*2) =P(I+NORD#2) 
CONTINUE 


DO 146 I=1,NORDD 
TMAX (I ) =AMAX1(TMAX(I) ,PRMT(I) ) 
TMIN(I)=AMIN1(TMIN(I),PRMT(I)) 
WRITE(3,200) (PRMT(J) ,J=1,NORDD) 
FORMAT (1X, 15F 12.4) 


UPDATE THE GAIN MATRIX P22(k) 
DO 120 I=1,NORDD 

DOMM22 J=1,NORDD 

P22 (1 J) =P22(1 J )-S(C1) #8 GJ) AM} 


MOVE U(k+1) TO U(k) 
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DO 124 L=1,NORD 
UP(LY1)=UN(L,1) 


CONTINUE 
CONTINUE 
WRITE(3,200) (TMAX(J) ,J=1,NORDD) 
WRITE(3,200) (TMIN(J) ,J=1,NORDD) 


WRITE(3, 300) 
FORMAT (1X, "ESTIMATE OF C VECTOR IS FOLLOWED BY 


So THAT OF A!) 


WRITE(3,301) 
FORMAT (1X,'IN FORMAT F12.4') 
WRITE(3,302) 
FORMAT(' ') 


WRITE (3,305) 
FORMAT(1X,'THESE ARE OUTPUT OF RML') 


RETURN 


END 


SUBROUTINE RAIBLE(A,N,ITEST) 


THIS PROGRAM TESTS THE STABILITY OF THE POLYNOMIAL 
F(Z) IN THE FORM 


N N=] 

Th Gh Ca NG ay A EP creiteced ba enn) 
A: AN ARRAY CONTAINING THE COEFFICIENTS OF F(Z) 
N: THE DEGREE OF THE POLYNOMIAL F(Z) 


LEST Ow hc) LoNoTADLE 


1 F(Z) IS UNSTABLE OR SIGULAR CASE 


REAL A(N),AA(6),BB(5) 
NN=N+1 


DO. 10: 1=2,NN 
AA(1I)=A(1I-1) 


TH GAWOIIOT Zt AST % hte 
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30 


40 


70 


ANCA) =a..0 

DO 70 K=1,N 
AK=AA (NN) /AA(1) 
NN=NN- 1 


DO 30 I=1,NN 
BB(I)=-AK#AA (NN-I+2) 


DO 40 I=1,NN 
AA(I)=AA(I)+BB(I) 


TEA AG eG T..0..0), GOTO: 70 
ITEST=1 

RETURN 

CONTINUE 

ITEST=0 


RETURN 
END 
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